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