Pruebas Pyomo

July 9, 2022

1 Pruebas con Pyomo

En este jupyter vamos a ir probando algunas cosas en pyomo, sobre como ir seteando variables y
como las va interpretando el paquete

2 cargado de los módulos

[1]: import [Link] as plt
import numpy as np
from [Link] import *
from [Link] import *

Vamos a hacer unos ejemplos sueltos para probar cosas

2.1 Ejemplos sueltos

En este caso vamos a declarar variables en un modelo para ver los valores y los reportes que se
obtienen
[2]: model_ejem = ConcreteModel()
model_ejem.A = Set(initialize = [1,2,3])
model_ejem.y = Var(within=model_ejem.A)
model_ejem.r = Var(domain=Reals)
model_ejem.w = Var(within=Boolean)

2.1.1 Otra forma de definir domain, via una función

[3]: def s_domain(model, i):

return RangeSet(i,i+1)
model_ejem.s = Var(model_ejem.A, domain=s_domain)

2.1.2 Imprimiendo el modelo

[4]: model_ejem.pprint()

1 Set Declarations
A : Size=1, Index=None, Ordered=Insertion

1
 Key : Dimen : Domain : Size : Members
None : 1 : Any : 3 : {1, 2, 3}

4 Var Declarations
r : Size=1, Index=None
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : None : None : None : False : True : Reals
s : Size=3, Index=A
Key : Lower : Value : Upper : Fixed : Stale : Domain
1 : 1 : None : 2 : False : True : [1:2]
2 : 2 : None : 3 : False : True : [2:3]
3 : 3 : None : 4 : False : True : [3:4]
w : Size=1, Index=None
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : 0 : None : 1 : False : True : Boolean
y : Size=1, Index=None
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : 1 : None : 3 : False : True : A

5 Declarations: A y r w s
En este ejemplo la variable s es ina variable indexada, cuyas componentes se definen sobre los
enteros consecutivos.

2.1.3 Caso en el que usamos funciones para los boundas, y inicialización usando dic-
tionary

[5]: lower = {1:2.5, 2:4.5, 3:6.5}

upper = {1:3.5 ,2:5.5 ,3:7.5}
init = {1:3., 2:5., 3:7.}
def b_bound(model, i):
return (lower[i], upper[i]) # return a 2-tuple with lower and upper bounds
def b_init(model, i):
return init[i] # return the corresponding index initialization value
model_ejem.b = Var(model_ejem.A, bounds=b_bound, initialize=b_init)

[6]: model_ejem.pprint()

1 Set Declarations
A : Size=1, Index=None, Ordered=Insertion
Key : Dimen : Domain : Size : Members
None : 1 : Any : 3 : {1, 2, 3}

5 Var Declarations
b : Size=3, Index=A
Key : Lower : Value : Upper : Fixed : Stale : Domain
1 : 2.5 : 3.0 : 3.5 : False : False : Reals
2 : 4.5 : 5.0 : 5.5 : False : False : Reals

2
 3 : 6.5 : 7.0 : 7.5 : False : False : Reals
r : Size=1, Index=None
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : None : None : None : False : True : Reals
s : Size=3, Index=A
Key : Lower : Value : Upper : Fixed : Stale : Domain
1 : 1 : None : 2 : False : True : [1:2]
2 : 2 : None : 3 : False : True : [2:3]
3 : 3 : None : 4 : False : True : [3:4]
w : Size=1, Index=None
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : 0 : None : 1 : False : True : Boolean
y : Size=1, Index=None
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : 1 : None : 3 : False : True : A

6 Declarations: A y r w s b

2.1.4 Otra forma

[7]: def g(model, i):
return 3*i
model_ejem.m = Var(model_ejem.A, initialize=g)

[8]: print(value(model_ejem.m[:]))

[3, 6, 9]

[9]: model_ejem.pprint()

1 Set Declarations
A : Size=1, Index=None, Ordered=Insertion
Key : Dimen : Domain : Size : Members
None : 1 : Any : 3 : {1, 2, 3}

6 Var Declarations
b : Size=3, Index=A
Key : Lower : Value : Upper : Fixed : Stale : Domain
1 : 2.5 : 3.0 : 3.5 : False : False : Reals
2 : 4.5 : 5.0 : 5.5 : False : False : Reals
3 : 6.5 : 7.0 : 7.5 : False : False : Reals
m : Size=3, Index=A
Key : Lower : Value : Upper : Fixed : Stale : Domain
1 : None : 3 : None : False : False : Reals
2 : None : 6 : None : False : False : Reals
3 : None : 9 : None : False : False : Reals
r : Size=1, Index=None
Key : Lower : Value : Upper : Fixed : Stale : Domain

3
 None : None : None : None : False : True : Reals
s : Size=3, Index=A
Key : Lower : Value : Upper : Fixed : Stale : Domain
1 : 1 : None : 2 : False : True : [1:2]
2 : 2 : None : 3 : False : True : [2:3]
3 : 3 : None : 4 : False : True : [3:4]
w : Size=1, Index=None
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : 0 : None : 1 : False : True : Boolean
y : Size=1, Index=None
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : 1 : None : 3 : False : True : A

7 Declarations: A y r w s b m

2.2 Usando el ejemplo del Libro/Presentación

Se requiere optimizar el siguiente problema

min 𝑥3 (𝑡𝑓 )
s.t. 𝑥1̇ = 𝑥2
𝑥2̇ = 𝑥2 + 𝑢 (1)
𝑥3̇ = 𝑥21 + 𝑥22 + 0.005 ∗ 𝑢2
𝑥2 − 8(𝑡 − 0.5)2 + 0.5 ≤ 0
𝑥1 (0) = 0, 𝑥2 (0) = −1, 𝑥3 (0) = 0, 𝑡𝑓 = 1

2.2.1 Declarando el modelo

[10]: # The concrete model object is instantiated
m = ConcreteModel()

# The parameter and time set

[Link] = Param(initialize=1)
m.t = ContinuousSet(bounds=(0, [Link]))

# The control variable

def u_init(m,t):
if t<0.15:
return 7.
elif t<0.3:
return 0.
elif t<0.5:
return -1.0
elif t<0.9:
return 1.
else:
return 0.
m.u = Var(m.t, initialize=u_init)

4
 # The state variables
m.x1 = Var(m.t)
m.x2 = Var(m.t)
m.x3 = Var(m.t)
# Derivatives of the state variable
m.dx1 = DerivativeVar(m.x1, wrt=m.t)
m.dx2 = DerivativeVar(m.x2, wrt=m.t)
m.dx3 = DerivativeVar(m.x3) # Assume the only varible to respec to

[11]: def plotter(subplot, x, *y, **kwds):

[Link](subplot)
for i, _y in enumerate(y):
[Link](list(x),[value(_y[t]) for t in x], 'brgcmk'[i%6])
if [Link]('points', False):
[Link](list(x), [value(_y[t]) for t in x], 'o')
[Link]([Link]('title',''))
[Link](tuple(_y.name for _y in y))
[Link]([Link])

fig, axes = [Link](nrows=1, ncols=1, figsize=(8, 8));

plotter(111, m.t, m.u, title='Initial Control Variable', points=True)

5
[12]: # First ODE
def _x1dot(m,t):
if t== [Link](): # override to initial conditions
return [Link]
return m.dx1[t] == m.x2[t]
m.x1dotcon = Constraint(m.t, rule=_x1dot)

# Second ODE
def _x2dot(m,t):
if t== [Link](): # override to initial conditions
return [Link]

6
 return m.dx2[t] == -m.x2[t] + m.u[t]
m.x2dotcon = Constraint(m.t, rule=_x2dot)

# Third ODE
def _x3dot(m,t):
if t== [Link](): # override to initial conditions
return [Link]
return m.dx3[t] == m.x1[t]**2 + m.x2[t]**2 + 0.005*m.u[t]**2
m.x3dotcon = Constraint(m.t, rule=_x3dot)

# The initial conditions as contrains (here as ConstrainList)

def _init(m):
#. Pay special attention that is used "==" instead of "="
yield m.x1[0] == 0. # using directly the boundary value of the continuous␣
↪set t

yield m.x2[[Link]()] == -1 # using the fisrt accesor function

yield m.x3[0] == 0
m.init_conditions = ConstraintList(rule = _init)

[13]: # Objective function

[Link] = Objective(expr=m.x3[[Link]]) # I think it could be used rule also see␣
↪the book

# The inequality constrain

def _con(m,t):
return m.x2[t] - 8*(t-0.5)**2 + 0.5 <= 0
[Link] = Constraint(m.t, rule=_con)

[14]: [Link]()

1 Set Declarations
init_conditions_index : Size=1, Index=None, Ordered=Insertion
Key : Dimen : Domain : Size : Members
None : 1 : Any : 3 : {1, 2, 3}

1 RangeSet Declarations
t_domain : Dimen=1, Size=Inf, Bounds=(0, 1)
Key : Finite : Members
None : False : [0..1]

1 Param Declarations
tf : Size=1, Index=None, Domain=Any, Default=None, Mutable=False
Key : Value
None : 1

4 Var Declarations
u : Size=2, Index=t

7
 Key : Lower : Value : Upper : Fixed : Stale : Domain
0 : None : 7.0 : None : False : False : Reals
1 : None : 0.0 : None : False : False : Reals
x1 : Size=2, Index=t
Key : Lower : Value : Upper : Fixed : Stale : Domain
0 : None : None : None : False : True : Reals
1 : None : None : None : False : True : Reals
x2 : Size=2, Index=t
Key : Lower : Value : Upper : Fixed : Stale : Domain
0 : None : None : None : False : True : Reals
1 : None : None : None : False : True : Reals
x3 : Size=2, Index=t
Key : Lower : Value : Upper : Fixed : Stale : Domain
0 : None : None : None : False : True : Reals
1 : None : None : None : False : True : Reals

1 Objective Declarations
obj : Size=1, Index=None, Active=True
Key : Active : Sense : Expression
None : True : minimize : x3[1]

5 Constraint Declarations
con : Size=2, Index=t, Active=True
Key : Lower : Body : Upper : Active
0 : -Inf : x2[0] - 2.0 + 0.5 : 0.0 : True
1 : -Inf : x2[1] - 2.0 + 0.5 : 0.0 : True
init_conditions : Size=3, Index=init_conditions_index, Active=True
Key : Lower : Body : Upper : Active
1 : 0.0 : x1[0] : 0.0 : True
2 : -1.0 : x2[0] : -1.0 : True
3 : 0.0 : x3[0] : 0.0 : True
x1dotcon : Size=1, Index=t, Active=True
Key : Lower : Body : Upper : Active
1 : 0.0 : dx1[1] - x2[1] : 0.0 : True
x2dotcon : Size=1, Index=t, Active=True
Key : Lower : Body : Upper : Active
1 : 0.0 : dx2[1] - (- x2[1] + u[1]) : 0.0 : True
x3dotcon : Size=1, Index=t, Active=True
Key : Lower : Body : Upper :
Active
1 : 0.0 : dx3[1] - (x1[1]**2 + x2[1]**2 + 0.005*u[1]**2) : 0.0 :
True

1 ContinuousSet Declarations
t : Size=1, Index=None, Ordered=Sorted
Key : Dimen : Domain : Size : Members
None : 1 : [0..1] : 2 : {0, 1}

8
 3 DerivativeVar Declarations
dx1 : Size=2, Index=t
Key : Lower : Value : Upper : Fixed : Stale : Domain
0 : None : None : None : False : True : Reals
1 : None : None : None : False : True : Reals
dx2 : Size=2, Index=t
Key : Lower : Value : Upper : Fixed : Stale : Domain
0 : None : None : None : False : True : Reals
1 : None : None : None : False : True : Reals
dx3 : Size=2, Index=t
Key : Lower : Value : Upper : Fixed : Stale : Domain
0 : None : None : None : False : True : Reals
1 : None : None : None : False : True : Reals

17 Declarations: tf t_domain t u x1 x2 x3 dx1 dx2 dx3 x1dotcon x2dotcon x3dotcon

init_conditions_index init_conditions obj con

2.2.2 Soolución del Modelo una vez declarado

El único método de solución para ODE disponible en [Link] es bajo el enfoque de discretización
simulatanea, tambien conocido como transcripción directa.
Este enfoque discretiza el dominio continuo in el modelo, y aproxima las EDOs usando ecua-
ciones algebraicas definidas en los puntos de discretización. El resultado de esta transformación
discretizada es un modelo puramente algebraico y es pasado para ser resuelto por un “solver” de
programación no-lineal estandar.
Hay dos tipos de esquemas de discretización: * Diferencia finita * Colocación
Esos esquemas difieren en las ecuaciones algebráicas usadas para aproximar la derivada, pero usan
escencialmente la misma sintáxis. Internamente, pyomo transforma el problema en un modelo
algebráico, y es lo que se conoce como transformación de modelo.

Diferencia Finita Para aplicar el método de diferencia finita

[15]: # # Construct the discretizer object
# discretizer = TransformationFactory('dae.finite_difference')
# # The parameters of the discretizer
# discretizer.apply_to(m, # the model to be applied
# nfe=20, # number of finite elements
# wrt=m.t, # with respect to what variable
# scheme='BACKWARD') # the used scheme also 'CENTRAL',␣
↪'FORWARD'

Colocación Para aplicar el método de colocación ortogonal sobre elmentos finitos

[16]: # Construct the discretizer object
discretizer2 = TransformationFactory('[Link]')
# The parameters of the discretizer to be applied

9
 discretizer2.apply_to(m, # the model to be applied
nfe=7, # number of finite elements
ncp=6, # number of collocation points
scheme='LAGRANGE-RADAU') # the used scheme also␣
↪'LAGRANGE-LEGENDRE'

Reduciendo el numero de grados de liberta en los puntos de colocación Para reducir los
grados de libertad, en una variable, en este caso la de control, debe usarse el comando a continuación
[17]: # If you want TO REDUCE the number of freedom inside of the
# finite elment, in a particular variable, you can do
discretizer2.reduce_collocation_points(m, # the model being applied
var=m.u, # the variable to be reduced
ncp=1, # the new number of␣
↪collocation points

contset=m.t) # respecting to which␣

↪continuous set

[17]: <[Link] at 0x7faf481c7680>

[18]: #[Link]()

[19]: def plotter(subplot, x, *y, **kwds):

[Link](subplot)
for i, _y in enumerate(y):
[Link](list(x),[value(_y[t]) for t in x], 'brgcmk'[i%6])
if [Link]('points', False):
[Link](list(x), [value(_y[t]) for t in x], 'o')
[Link]([Link]('title',''))
[Link](tuple(_y.name for _y in y))
[Link]([Link])

fig, axes = [Link](nrows=1, ncols=1, figsize=(8, 8));

plotter(111, m.t, m.u, title='Initial Control Variable', points=True)

10
 2.2.3 Aplicando el solver
En este caso se va a usar el solver del puntos interiores
[20]: solver = SolverFactory('ipopt')
results = [Link](m, tee=True).write()

Ipopt 3.14.7:

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit [Link]

11
******************************************************************************

This is Ipopt version 3.14.7, running with linear solver MUMPS 5.2.1.

Number of nonzeros in equality constraint Jacobian…: 1459

Number of nonzeros in inequality constraint Jacobian.: 43
Number of nonzeros in Lagrangian Hessian…: 126

Total number of variables…: 297

variables with only lower bounds: 0
variables with lower and upper bounds: 0
variables with only upper bounds: 0
Total number of equality constraints…: 290
Total number of inequality constraints…: 43
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 43

iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls

0 0.0000000e+00 7.00e+00 2.71e-01 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 3.5766270e-02 9.10e-01 4.56e+00 -1.0 9.33e+00 - 1.78e-01 1.00e+00f 1
2 3.1480130e-01 2.72e-02 6.18e-01 -1.0 9.76e-01 - 6.89e-01 1.00e+00h 1
3 3.6879902e-01 5.64e-03 2.35e-02 -1.7 4.95e-01 - 9.49e-01 1.00e+00h 1
4 2.6591821e-01 3.22e-02 2.83e-08 -2.5 7.70e-01 - 1.00e+00 1.00e+00h 1
5 2.4034588e-01 2.85e-03 2.83e-08 -2.5 2.26e-01 - 1.00e+00 1.00e+00h 1
6 1.9151649e-01 1.07e-02 5.89e-04 -3.8 4.18e-01 - 1.00e+00 9.32e-01h 1
7 1.9194009e-01 2.71e-03 2.44e-02 -3.8 2.83e-01 - 1.00e+00 8.68e-01h 1
8 1.9348671e-01 4.45e-05 1.50e-09 -3.8 6.62e-02 - 1.00e+00 1.00e+00h 1
9 1.9248019e-01 3.99e-05 7.44e-03 -5.7 3.95e-02 - 1.00e+00 7.95e-01h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 1.9246497e-01 2.59e-06 1.84e-11 -5.7 1.13e-02 - 1.00e+00 1.00e+00h 1
11 1.9246119e-01 6.89e-08 3.49e-06 -8.6 2.57e-03 - 9.92e-01 9.96e-01h 1
12 1.9246121e-01 1.06e-10 2.51e-14 -8.6 1.04e-04 - 1.00e+00 1.00e+00h 1

Number of Iterations…: 12

(scaled) (unscaled)
Objective…: 1.9246120821985307e-01 1.9246120821985307e-01
Dual infeasibility…: 2.5063284780912909e-14 2.5063284780912909e-14
Constraint violation…: 1.0606591993589376e-10 1.0606591993589376e-10
Variable bound violation: 0.0000000000000000e+00 0.0000000000000000e+00
Complementarity…: 2.5272282074814942e-09 2.5272282074814942e-09
Overall NLP error…: 2.5272282074814942e-09 2.5272282074814942e-09

Number of objective function evaluations = 13

Number of objective gradient evaluations = 13
Number of equality constraint evaluations = 13

12
 Number of inequality constraint evaluations = 13
Number of equality constraint Jacobian evaluations = 13
Number of inequality constraint Jacobian evaluations = 13
Number of Lagrangian Hessian evaluations = 12
Total seconds in IPOPT = 0.229

EXIT: Optimal Solution Found.

# ==========================================================
# = Solver Results =
# ==========================================================
# ----------------------------------------------------------
# Problem Information
# ----------------------------------------------------------
Problem:
- Lower bound: -inf
Upper bound: inf
Number of objectives: 1
Number of constraints: 333
Number of variables: 297
Sense: unknown
# ----------------------------------------------------------
# Solver Information
# ----------------------------------------------------------
Solver:
- Status: ok
Message: Ipopt 3.14.7\x3a Optimal Solution Found
Termination condition: optimal
Id: 0
Error rc: 0
Time: 0.27477407455444336
# ----------------------------------------------------------
# Solution Information
# ----------------------------------------------------------
Solution:
- number of solutions: 0
number of solutions displayed: 0

[21]: def plotter(subplot, x, *y, **kwds):

[Link](subplot)
for i, _y in enumerate(y):
[Link](list(x),[value(_y[t]) for t in x], 'brgcmk'[i%6])
if [Link]('points', False):
[Link](list(x), [value(_y[t]) for t in x], 'o')
[Link]([Link]('title',''))
[Link](tuple(_y.name for _y in y))
[Link]([Link])

13
 fig, axes = [Link](nrows=1, ncols=2, figsize=(14, 8));
plotter(121, m.t, m.x1, m.x2, title='Differential Variables')
plotter(122, m.t, m.u, title='Control Variable', points=True)

[22]: print( value(m.u[:]))

[7.0, 6.889205315426058, 6.889205315426058, 6.889205315426058,

6.889205315426058, 6.889205315426058, 6.889205315426058, -1.5587072215166708,
-1.5587072215166708, -1.5587072215166708, -1.5587072215166708,
-1.5587072215166708, -1.5587072215166708, -2.668390704569762,
-2.668390704569762, -2.668390704569762, -2.668390704569762, -2.668390704569762,
-2.668390704569762, -0.5306411266010939, -0.5306411266010939,
-0.5306411266010939, -0.5306411266010939, -0.5306411266010939,
-0.5306411266010939, 1.8959807787484528, 1.8959807787484528, 1.8959807787484528,
1.8959807787484528, 1.8959807787484528, 1.8959807787484528, 1.3388461735804136,
1.3388461735804136, 1.3388461735804136, 1.3388461735804136, 1.3388461735804136,
1.3388461735804136, -0.03699436632100241, -0.03699436632100241,
-0.03699436632100241, -0.03699436632100241, -0.03699436632100241,
-0.03699436632100241]

[ ]:

14