1# ___________________________________________________________________________ 2# 3# Pyomo: Python Optimization Modeling Objects 4# Copyright 2017 National Technology and Engineering Solutions of Sandia, LLC 5# Under the terms of Contract DE-NA0003525 with National Technology and 6# Engineering Solutions of Sandia, LLC, the U.S. Government retains certain 7# rights in this software. 8# This software is distributed under the 3-clause BSD License. 9# ___________________________________________________________________________ 10 11import pyomo.environ as pyo 12import pyomo.dae as dae 13 14 15def make_gas_expansion_model(N=2): 16 """ 17 This is the simplest model I could think of that has a 18 subsystem with a non-trivial block triangularization. 19 Something like a gas (somehow) undergoing a series 20 of isentropic expansions. 21 """ 22 m = pyo.ConcreteModel() 23 m.streams = pyo.Set(initialize=range(N+1)) 24 m.rho = pyo.Var(m.streams, initialize=1) 25 m.P = pyo.Var(m.streams, initialize=1) 26 m.F = pyo.Var(m.streams, initialize=1) 27 m.T = pyo.Var(m.streams, initialize=1) 28 29 m.R = pyo.Param(initialize=8.31) 30 m.Q = pyo.Param(m.streams, initialize=1) 31 m.gamma = pyo.Param(initialize=1.4*m.R.value) 32 33 def mbal(m, i): 34 if i == 0: 35 return pyo.Constraint.Skip 36 else: 37 return m.rho[i-1]*m.F[i-1] - m.rho[i]*m.F[i] == 0 38 m.mbal = pyo.Constraint(m.streams, rule=mbal) 39 40 def ebal(m, i): 41 if i == 0: 42 return pyo.Constraint.Skip 43 else: 44 return ( 45 m.rho[i-1]*m.F[i-1]*m.T[i-1] + 46 m.Q[i] - 47 m.rho[i]*m.F[i]*m.T[i] == 0 48 ) 49 m.ebal = pyo.Constraint(m.streams, rule=ebal) 50 51 def expansion(m, i): 52 if i == 0: 53 return pyo.Constraint.Skip 54 else: 55 return m.P[i]/m.P[i-1] - (m.rho[i]/m.rho[i-1])**m.gamma == 0 56 m.expansion = pyo.Constraint(m.streams, rule=expansion) 57 58 def ideal_gas(m, i): 59 return m.P[i] - m.rho[i]*m.R*m.T[i] == 0 60 m.ideal_gas = pyo.Constraint(m.streams, rule=ideal_gas) 61 62 return m 63 64 65def make_dynamic_model(**disc_args): 66 # Level control model 67 m = pyo.ConcreteModel() 68 m.time = dae.ContinuousSet(initialize=[0.0, 10.0]) 69 m.height = pyo.Var(m.time, initialize=1.0) 70 m.flow_in = pyo.Var(m.time, initialize=1.0) 71 m.flow_out = pyo.Var(m.time, initialize=0.5) 72 m.dhdt = dae.DerivativeVar(m.height, wrt=m.time, initialize=0.0) 73 74 m.area = pyo.Param(initialize=1.0) 75 m.flow_const = pyo.Param(initialize=0.5) 76 77 def diff_eqn_rule(m, t): 78 return m.area*m.dhdt[t] - (m.flow_in[t] - m.flow_out[t]) == 0 79 m.diff_eqn = pyo.Constraint(m.time, rule=diff_eqn_rule) 80 81 def flow_out_rule(m, t): 82 return m.flow_out[t] - (m.flow_const*pyo.sqrt(m.height[t])) == 0 83 m.flow_out_eqn = pyo.Constraint(m.time, rule=flow_out_rule) 84 85 default_disc_args = { 86 "wrt": m.time, 87 "nfe": 5, 88 "scheme": "BACKWARD", 89 } 90 default_disc_args.update(disc_args) 91 92 discretizer = pyo.TransformationFactory("dae.finite_difference") 93 discretizer.apply_to(m, **default_disc_args) 94 95 return m 96 97 98def make_degenerate_solid_phase_model(): 99 """ 100 From the solid phase thermo package of a moving bed chemical looping 101 combustion reactor. This example was first presented in [1] 102 103 [1] Parker, R. Nonlinear programming strategies for dynamic models of 104 chemical looping combustion reactors. Pres. AIChE Annual Meeting, 2020. 105 106 """ 107 m = pyo.ConcreteModel() 108 m.components = pyo.Set(initialize=[1,2,3]) 109 m.x = pyo.Var(m.components, initialize=1/3) 110 m.flow_comp = pyo.Var(m.components, initialize=10) 111 m.flow = pyo.Var(initialize=30) 112 m.rho = pyo.Var(initialize=1) 113 114 # These are rough approximations of the relevant equations, with the same 115 # incidence. 116 m.sum_eqn = pyo.Constraint(expr=sum(m.x[j] for j in m.components) - 1 == 0) 117 m.holdup_eqn = pyo.Constraint(m.components, expr={ 118 j: m.x[j]*m.rho - 1 == 0 for j in m.components 119 }) 120 m.density_eqn = pyo.Constraint(expr= 121 1/m.rho - sum(1/m.x[j] for j in m.components) == 0 122 ) 123 m.flow_eqn = pyo.Constraint(m.components, expr={ 124 j: m.x[j]*m.flow - m.flow_comp[j] == 0 for j in m.components 125 }) 126 127 return m 128