Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Apply black to examples directory py files #2720

Merged
merged 2 commits into from
Feb 10, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
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