Skip to content

Commit

Permalink
Apply black to the examples directory py files
Browse files Browse the repository at this point in the history
  • Loading branch information
mrmundt committed Feb 7, 2023
1 parent 9ee4235 commit 3fa9f89
Show file tree
Hide file tree
Showing 295 changed files with 9,069 additions and 4,472 deletions.
58 changes: 35 additions & 23 deletions examples/dae/Heat_Conduction.py
Original file line number Diff line number Diff line change
Expand Up @@ -15,37 +15,49 @@
from pyomo.dae import *

m = ConcreteModel()
m.time = ContinuousSet(bounds=(0,1))
m.x = ContinuousSet(bounds=(0,10))
m.y = ContinuousSet(bounds=(0,5))
m.T = Var(m.x,m.y,m.time)
m.u = Var(m.x,m.y,m.time)
m.time = ContinuousSet(bounds=(0, 1))
m.x = ContinuousSet(bounds=(0, 10))
m.y = ContinuousSet(bounds=(0, 5))
m.T = Var(m.x, m.y, m.time)
m.u = Var(m.x, m.y, m.time)
m.T0 = Param(initialize=5)
m.TD = Param(m.x,m.y,initialize=25)
m.TD = Param(m.x, m.y, initialize=25)
m.Ux0 = Param(initialize=10)
m.Uy5 = Param(initialize=15)

m.dTdx = DerivativeVar(m.T,wrt=m.x)
m.d2Tdx2 = DerivativeVar(m.T,wrt=(m.x,m.x))
m.dTdy = DerivativeVar(m.T,wrt=m.y)
m.d2Tdy2 = DerivativeVar(m.T,wrt=(m.y,m.y))
m.dTdt = DerivativeVar(m.T,wrt=m.time)
m.dTdx = DerivativeVar(m.T, wrt=m.x)
m.d2Tdx2 = DerivativeVar(m.T, wrt=(m.x, m.x))
m.dTdy = DerivativeVar(m.T, wrt=m.y)
m.d2Tdy2 = DerivativeVar(m.T, wrt=(m.y, m.y))
m.dTdt = DerivativeVar(m.T, wrt=m.time)

def _heateq(m,i,j,k):
return m.d2Tdx2[i,j,k] + m.d2Tdy2[i,j,k] + m.u[i,j,k] == m.dTdt[i,j,k]
m.heateq = Constraint(m.x,m.y,m.time,rule=_heateq)

def _initT(m,i,j):
return m.T[i,j,0] == m.T0
m.initT = Constraint(m.x,m.y,rule=_initT)
def _heateq(m, i, j, k):
return m.d2Tdx2[i, j, k] + m.d2Tdy2[i, j, k] + m.u[i, j, k] == m.dTdt[i, j, k]

def _xbound(m,j,k):
return m.dTdx[0,j,k] == m.Ux0
m.xbound = Constraint(m.y,m.time,rule=_xbound)

def _ybound(m,i,k):
return m.dTdy[i,5,k] == m.Uy5
m.ybound = Constraint(m.x,m.time,rule=_ybound)
m.heateq = Constraint(m.x, m.y, m.time, rule=_heateq)


def _initT(m, i, j):
return m.T[i, j, 0] == m.T0


m.initT = Constraint(m.x, m.y, rule=_initT)


def _xbound(m, j, k):
return m.dTdx[0, j, k] == m.Ux0


m.xbound = Constraint(m.y, m.time, rule=_xbound)


def _ybound(m, i, k):
return m.dTdy[i, 5, k] == m.Uy5


m.ybound = Constraint(m.x, m.time, rule=_ybound)

# def _intExp(m,i,j):
# return m.T[i,j,1] - m.TD[i,j]
Expand Down
46 changes: 27 additions & 19 deletions examples/dae/Optimal_Control.py
Original file line number Diff line number Diff line change
Expand Up @@ -11,42 +11,50 @@

# Sample Problem 1 (Ex 1 from Dynopt Guide)
#
# min X2(tf)
# s.t. X1_dot = u X1(0) = 1
# X2_dot = X1^2 + u^2 X2(0) = 0
# tf = 1
# min X2(tf)
# s.t. X1_dot = u X1(0) = 1
# X2_dot = X1^2 + u^2 X2(0) = 0
# tf = 1

from pyomo.environ import *
from pyomo.dae import *

m = ConcreteModel()

m.t = ContinuousSet(bounds=(0,1))
m.t = ContinuousSet(bounds=(0, 1))

m.x1 = Var(m.t, bounds=(0,1))
m.x2 = Var(m.t, bounds=(0,1))
m.x1 = Var(m.t, bounds=(0, 1))
m.x2 = Var(m.t, bounds=(0, 1))
m.u = Var(m.t, initialize=0)

m.x1dot = DerivativeVar(m.x1)
m.x2dot = DerivativeVar(m.x2)

m.obj = Objective(expr=m.x2[1])

def _x1dot(M,i):
if i == 0:
return Constraint.Skip
return M.x1dot[i] == M.u[i]

def _x1dot(M, i):
if i == 0:
return Constraint.Skip
return M.x1dot[i] == M.u[i]


m.x1dotcon = Constraint(m.t, rule=_x1dot)

def _x2dot(M,i):
if i == 0:
return Constraint.Skip
return M.x2dot[i] == M.x1[i]**2 + M.u[i]**2

def _x2dot(M, i):
if i == 0:
return Constraint.Skip
return M.x2dot[i] == M.x1[i] ** 2 + M.u[i] ** 2


m.x2dotcon = Constraint(m.t, rule=_x2dot)


def _init(M):
yield M.x1[0] == 1
yield M.x2[0] == 0
yield ConstraintList.End
m.init_conditions = ConstraintList(rule=_init)
yield M.x1[0] == 1
yield M.x2[0] == 0
yield ConstraintList.End


m.init_conditions = ConstraintList(rule=_init)
67 changes: 40 additions & 27 deletions examples/dae/PDE_example.py
Original file line number Diff line number Diff line change
Expand Up @@ -16,33 +16,45 @@

m = ConcreteModel()
m.pi = Param(initialize=3.1416)
m.t = ContinuousSet(bounds=(0,2))
m.x = ContinuousSet(bounds=(0,1))
m.u = Var(m.x,m.t)
m.t = ContinuousSet(bounds=(0, 2))
m.x = ContinuousSet(bounds=(0, 1))
m.u = Var(m.x, m.t)

m.dudx = DerivativeVar(m.u,wrt=m.x)
m.dudx2 = DerivativeVar(m.u,wrt=(m.x,m.x))
m.dudt = DerivativeVar(m.u,wrt=m.t)
m.dudx = DerivativeVar(m.u, wrt=m.x)
m.dudx2 = DerivativeVar(m.u, wrt=(m.x, m.x))
m.dudt = DerivativeVar(m.u, wrt=m.t)

def _pde(m,i,j):
if i == 0 or i == 1 or j == 0 :

def _pde(m, i, j):
if i == 0 or i == 1 or j == 0:
return Constraint.Skip
return m.pi**2*m.dudt[i,j] == m.dudx2[i,j]
m.pde = Constraint(m.x,m.t,rule=_pde)
return m.pi**2 * m.dudt[i, j] == m.dudx2[i, j]


m.pde = Constraint(m.x, m.t, rule=_pde)


def _initcon(m,i):
def _initcon(m, i):
if i == 0 or i == 1:
return Constraint.Skip
return m.u[i,0] == sin(m.pi*i)
m.initcon = Constraint(m.x,rule=_initcon)
return m.u[i, 0] == sin(m.pi * i)


m.initcon = Constraint(m.x, rule=_initcon)


def _lowerbound(m,j):
return m.u[0,j] == 0
m.lowerbound = Constraint(m.t,rule=_lowerbound)
def _lowerbound(m, j):
return m.u[0, j] == 0

def _upperbound(m,j):
return m.pi*exp(-j)+m.dudx[1,j] == 0
m.upperbound = Constraint(m.t,rule=_upperbound)

m.lowerbound = Constraint(m.t, rule=_lowerbound)


def _upperbound(m, j):
return m.pi * exp(-j) + m.dudx[1, j] == 0


m.upperbound = Constraint(m.t, rule=_upperbound)

m.obj = Objective(expr=1)

Expand All @@ -54,27 +66,27 @@ def _upperbound(m,j):
# Discretize using Finite Difference and Collocation
discretizer = TransformationFactory('dae.finite_difference')
discretizer2 = TransformationFactory('dae.collocation')
discretizer.apply_to(m,nfe=25,wrt=m.x,scheme='BACKWARD')
discretizer2.apply_to(m,nfe=20,ncp=3,wrt=m.t)
discretizer.apply_to(m, nfe=25, wrt=m.x, scheme='BACKWARD')
discretizer2.apply_to(m, nfe=20, ncp=3, wrt=m.t)

# Discretize using Finite Difference Method
# discretizer = TransformationFactory('dae.finite_difference')
# discretizer.apply_to(m,nfe=25,wrt=m.x,scheme='BACKWARD')
# discretizer.apply_to(m,nfe=20,wrt=m.t,scheme='BACKWARD')

solver=SolverFactory('ipopt')
results = solver.solve(m,tee=True)
solver = SolverFactory('ipopt')
results = solver.solve(m, tee=True)

x = []
t = []
u = []

for i in sorted(m.x):
temp=[]
temp = []
tempx = []
for j in sorted(m.t):
tempx.append(i)
temp.append(value(m.u[i,j]))
temp.append(value(m.u[i, j]))
x.append(tempx)
t.append(sorted(m.t))
u.append(temp)
Expand All @@ -83,9 +95,10 @@ def _upperbound(m,j):
import numpy
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.axes3d import Axes3D

fig = plt.figure()
ax = fig.add_subplot(1,1,1,projection='3d')
ax = fig.add_subplot(1, 1, 1, projection='3d')
ax.set_xlabel('Distance x')
ax.set_ylabel('Time t')
p = ax.plot_wireframe(x,t,u,rstride=1,cstride=1)
p = ax.plot_wireframe(x, t, u, rstride=1, cstride=1)
fig.show()
58 changes: 35 additions & 23 deletions examples/dae/Parameter_Estimation.py
Original file line number Diff line number Diff line change
Expand Up @@ -12,52 +12,64 @@
# Sample Problem 2: Parameter Estimation
# (Ex 5 from Dynopt Guide)
#
# min sum((X1(ti)-X1_meas(ti))^2)
# s.t. X1_dot = X2 X1(0) = p1
# X2_dot = 1-2*X2-X1 X2(0) = p2
# -1.5 <= p1,p2 <= 1.5
# tf = 6
# min sum((X1(ti)-X1_meas(ti))^2)
# s.t. X1_dot = X2 X1(0) = p1
# X2_dot = 1-2*X2-X1 X2(0) = p2
# -1.5 <= p1,p2 <= 1.5
# tf = 6
#

from pyomo.environ import *
from pyomo.dae import *

model = AbstractModel()
model.t = ContinuousSet()
model.MEAS_t = Set(within=model.t) # Measurement times, must be subset of t
model.t = ContinuousSet()
model.MEAS_t = Set(within=model.t) # Measurement times, must be subset of t
model.x1_meas = Param(model.MEAS_t)

model.x1 = Var(model.t)
model.x2 = Var(model.t)

model.p1 = Var(bounds=(-1.5,1.5))
model.p2 = Var(bounds=(-1.5,1.5))
model.p1 = Var(bounds=(-1.5, 1.5))
model.p2 = Var(bounds=(-1.5, 1.5))

model.x1dot = DerivativeVar(model.x1,wrt=model.t)
model.x1dot = DerivativeVar(model.x1, wrt=model.t)
model.x2dot = DerivativeVar(model.x2)


def _init_conditions(model):
yield model.x1[0] == model.p1
yield model.x2[0] == model.p2
yield model.x1[0] == model.p1
yield model.x2[0] == model.p2


model.init_conditions = ConstraintList(rule=_init_conditions)

# Alternate way to declare initial conditions
#def _initx1(model):
# return model.x1[0] == model.p1
#model.initx1 = Constraint(rule=_initx1)
# def _initx1(model):
# return model.x1[0] == model.p1
# model.initx1 = Constraint(rule=_initx1)

# def _initx2(model):
# return model.x2[0] == model.p2
# model.initx2 = Constraint(rule=_initx2)


def _x1dot(model, i):
return model.x1dot[i] == model.x2[i]

#def _initx2(model):
# return model.x2[0] == model.p2
#model.initx2 = Constraint(rule=_initx2)

def _x1dot(model,i):
return model.x1dot[i] == model.x2[i]
model.x1dotcon = Constraint(model.t, rule=_x1dot)

def _x2dot(model,i):
return model.x2dot[i] == 1-2*model.x2[i]-model.x1[i]

def _x2dot(model, i):
return model.x2dot[i] == 1 - 2 * model.x2[i] - model.x1[i]


model.x2dotcon = Constraint(model.t, rule=_x2dot)


def _obj(model):
return sum((model.x1[i]-model.x1_meas[i])**2 for i in model.MEAS_t)
return sum((model.x1[i] - model.x1_meas[i]) ** 2 for i in model.MEAS_t)


model.obj = Objective(rule=_obj)
Loading

0 comments on commit 3fa9f89

Please sign in to comment.