0
# Differential Algebraic Equations
1

2
Components for modeling dynamic systems with differential and algebraic equations. DAE capabilities enable optimization of systems with time-dependent behavior, continuous processes, and dynamic constraints.
3

4
## Capabilities
5

6
### Continuous Set Components
7

8
Components for defining continuous domains and time discretization in dynamic optimization problems.
9

10
```python { .api }
11
class ContinuousSet:
12
    """
13
    Continuous set component for DAE modeling.
14
    
15
    Represents a continuous domain (typically time) that can be
16
    discretized for numerical solution of differential equations.
17
    """
18
    def __init__(self, *args, **kwargs): ...
19
    
20
    def get_finite_elements(self):
21
        """
22
        Get finite element discretization points.
23
        
24
        Returns:
25
            list: Discretization points
26
        """
27
    
28
    def get_discretization_info(self):
29
        """
30
        Get discretization scheme information.
31
        
32
        Returns:
33
            dict: Discretization metadata
34
        """
35
```
36

37
### Derivative Variable Components
38

39
Components for defining derivative variables and their relationships to state variables.
40

41
```python { .api }
42
class DerivativeVar:
43
    """
44
    Derivative variable component for differential equations.
45
    
46
    Represents the derivative of a state variable with respect
47
    to a continuous set (typically time).
48
    """
49
    def __init__(self, *args, **kwargs): ...
50
    
51
    def get_state_var(self):
52
        """
53
        Get the state variable this derivative represents.
54
        
55
        Returns:
56
            Var: Associated state variable
57
        """
58
    
59
    def get_continuousset(self):
60
        """
61
        Get the continuous set for differentiation.
62
        
63
        Returns:
64
            ContinuousSet: Continuous domain
65
        """
66
```
67

68
### Integral Components
69

70
Components for defining integral expressions and constraints in dynamic systems.
71

72
```python { .api }
73
class Integral:
74
    """
75
    Integral component for dynamic optimization.
76
    
77
    Represents integral expressions over continuous domains,
78
    commonly used in objective functions and constraints.
79
    """
80
    def __init__(self, *args, **kwargs): ...
81
    
82
    def get_integrand(self):
83
        """
84
        Get the integrand expression.
85
        
86
        Returns:
87
            Expression: Function being integrated
88
        """
89
    
90
    def get_continuousset(self):
91
        """
92
        Get the integration domain.
93
        
94
        Returns:
95
            ContinuousSet: Domain of integration
96
        """
97
```
98

99
### Simulation Interface
100

101
Interface for simulating DAE systems and initializing optimization models.
102

103
```python { .api }
104
class Simulator:
105
    """
106
    DAE simulation interface for model initialization and analysis.
107
    
108
    Provides capabilities for simulating differential-algebraic
109
    equation systems to generate initial conditions and validate models.
110
    """
111
    def __init__(self, model, **kwargs): ...
112
    
113
    def simulate(self, numpoints=None, integrator=None):
114
        """
115
        Simulate the DAE system.
116
        
117
        Args:
118
            numpoints (int, optional): Number of simulation points
119
            integrator (str, optional): Integration method
120
            
121
        Returns:
122
            dict: Simulation results
123
        """
124
    
125
    def initialize_model(self):
126
        """
127
        Initialize model with simulation results.
128
        """
129
    
130
    def get_profile(self, variable):
131
        """
132
        Get variable profile from simulation.
133
        
134
        Args:
135
            variable: Variable to extract profile for
136
            
137
        Returns:
138
            dict: Time-indexed variable values
139
        """
140
```
141

142
### DAE Error Handling
143

144
Exception class for DAE-specific errors and modeling issues.
145

146
```python { .api }
147
class DAE_Error(Exception):
148
    """DAE-specific error class for differential equation modeling issues."""
149
    def __init__(self, message): ...
150
```
151

152
## Usage Examples
153

154
### Basic DAE Model
155

156
```python
157
from pyomo.environ import *
158
from pyomo.dae import ContinuousSet, DerivativeVar
159

160
model = ConcreteModel()
161

162
# Define time domain
163
model.time = ContinuousSet(bounds=(0, 10))
164

165
# State variables
166
model.x1 = Var(model.time, bounds=(0, None))  # Concentration A
167
model.x2 = Var(model.time, bounds=(0, None))  # Concentration B
168

169
# Derivative variables
170
model.dx1dt = DerivativeVar(model.x1, wrt=model.time)
171
model.dx2dt = DerivativeVar(model.x2, wrt=model.time)
172

173
# Parameters
174
model.k1 = Param(initialize=2.5)  # Reaction rate constant 1
175
model.k2 = Param(initialize=1.0)  # Reaction rate constant 2
176

177
# Initial conditions
178
model.x1[0].fix(1.0)  # Initial concentration A
179
model.x2[0].fix(0.0)  # Initial concentration B
180

181
# Differential equations: A -> B -> C
182
def ode1_rule(model, t):
183
    return model.dx1dt[t] == -model.k1 * model.x1[t]
184

185
def ode2_rule(model, t):
186
    return model.dx2dt[t] == model.k1 * model.x1[t] - model.k2 * model.x2[t]
187

188
model.ode1 = Constraint(model.time, rule=ode1_rule)
189
model.ode2 = Constraint(model.time, rule=ode2_rule)
190

191
# Apply discretization
192
discretizer = TransformationFactory('dae.finite_difference')
193
discretizer.apply_to(model, nfe=20, wrt=model.time)
194
```
195

196
### Dynamic Optimization with Control
197

198
```python
199
from pyomo.environ import *
200
from pyomo.dae import ContinuousSet, DerivativeVar, Integral
201

202
model = ConcreteModel()
203

204
# Time domain
205
model.time = ContinuousSet(bounds=(0, 1))
206

207
# State variables
208
model.x1 = Var(model.time, bounds=(-10, 10))  # Position
209
model.x2 = Var(model.time, bounds=(-10, 10))  # Velocity
210

211
# Control variable
212
model.u = Var(model.time, bounds=(-1, 1))  # Control input
213

214
# Derivatives
215
model.dx1dt = DerivativeVar(model.x1, wrt=model.time)
216
model.dx2dt = DerivativeVar(model.x2, wrt=model.time)
217

218
# Initial conditions
219
model.x1[0].fix(0)  # Start at origin
220
model.x2[0].fix(0)  # Start at rest
221

222
# Terminal conditions
223
model.x1[1].fix(1)  # End at position 1
224
model.x2[1].fix(0)  # End at rest
225

226
# System dynamics
227
def position_ode(model, t):
228
    return model.dx1dt[t] == model.x2[t]
229

230
def velocity_ode(model, t):
231
    return model.dx2dt[t] == model.u[t]
232

233
model.position_ode = Constraint(model.time, rule=position_ode)
234
model.velocity_ode = Constraint(model.time, rule=velocity_ode)
235

236
# Objective: minimize control effort
237
model.obj_integrand = Var(model.time)
238
model.integrand_def = Constraint(
239
    model.time,
240
    rule=lambda model, t: model.obj_integrand[t] == model.u[t]**2
241
)
242

243
model.objective_integral = Integral(
244
    model.time,
245
    wrt=model.time,
246
    rule=lambda model, t: model.obj_integrand[t]
247
)
248

249
model.obj = Objective(expr=model.objective_integral, sense=minimize)
250

251
# Apply discretization
252
discretizer = TransformationFactory('dae.collocation')
253
discretizer.apply_to(model, nfe=10, ncp=3, wrt=model.time)
254
```
255

256
### Batch Reactor Optimization
257

258
```python
259
from pyomo.environ import *
260
from pyomo.dae import ContinuousSet, DerivativeVar, Integral
261

262
model = ConcreteModel()
263

264
# Time domain for batch process
265
model.time = ContinuousSet(bounds=(0, 2))  # 2 hours
266

267
# State variables
268
model.CA = Var(model.time, bounds=(0, 10))    # Concentration of A
269
model.CB = Var(model.time, bounds=(0, 10))    # Concentration of B
270
model.CC = Var(model.time, bounds=(0, 10))    # Concentration of C
271
model.T = Var(model.time, bounds=(300, 400))  # Temperature
272

273
# Control variable
274
model.Q = Var(model.time, bounds=(-50, 50))   # Heat input
275

276
# Derivatives
277
model.dCAdt = DerivativeVar(model.CA, wrt=model.time)
278
model.dCBdt = DerivativeVar(model.CB, wrt=model.time)
279
model.dCCdt = DerivativeVar(model.CC, wrt=model.time)
280
model.dTdt = DerivativeVar(model.T, wrt=model.time)
281

282
# Initial conditions
283
model.CA[0].fix(2.0)    # Initial concentration A
284
model.CB[0].fix(0.0)    # Initial concentration B  
285
model.CC[0].fix(0.0)    # Initial concentration C
286
model.T[0].fix(350)     # Initial temperature
287

288
# Parameters
289
model.k0 = Param(initialize=1e6)   # Pre-exponential factor
290
model.E = Param(initialize=8000)   # Activation energy
291
model.R = Param(initialize=8.314)  # Gas constant
292
model.rho = Param(initialize=1000) # Density
293
model.Cp = Param(initialize=1.0)   # Heat capacity
294

295
# Reaction rate expressions
296
def k_expr(model, t):
297
    return model.k0 * exp(-model.E / (model.R * model.T[t]))
298

299
def reaction_rate(model, t):
300
    return k_expr(model, t) * model.CA[t]
301

302
# Mass balance ODEs: A -> B -> C
303
def mass_balance_A(model, t):
304
    return model.dCAdt[t] == -reaction_rate(model, t)
305

306
def mass_balance_B(model, t):
307
    return model.dCBdt[t] == reaction_rate(model, t) - 0.5 * reaction_rate(model, t)
308

309
def mass_balance_C(model, t):
310
    return model.dCCdt[t] == 0.5 * reaction_rate(model, t)
311

312
# Energy balance
313
def energy_balance(model, t):
314
    return model.rho * model.Cp * model.dTdt[t] == \
315
           -10000 * reaction_rate(model, t) + model.Q[t]
316

317
model.mass_A = Constraint(model.time, rule=mass_balance_A)
318
model.mass_B = Constraint(model.time, rule=mass_balance_B)
319
model.mass_C = Constraint(model.time, rule=mass_balance_C)
320
model.energy = Constraint(model.time, rule=energy_balance)
321

322
# Objective: maximize final concentration of B
323
model.obj = Objective(expr=model.CB[2], sense=maximize)
324

325
# Apply discretization
326
discretizer = TransformationFactory('dae.finite_difference')
327
discretizer.apply_to(model, nfe=40, wrt=model.time, scheme='BACKWARD')
328
```
329

330
### DAE Model Simulation and Initialization
331

332
```python
333
from pyomo.environ import *
334
from pyomo.dae import ContinuousSet, DerivativeVar, Simulator
335

336
# Create DAE model
337
model = ConcreteModel()
338
model.time = ContinuousSet(bounds=(0, 5))
339
model.x = Var(model.time)
340
model.dxdt = DerivativeVar(model.x, wrt=model.time)
341

342
# Model equations
343
model.x[0].fix(1.0)
344
model.ode = Constraint(
345
    model.time,
346
    rule=lambda m, t: m.dxdt[t] == -0.5 * m.x[t]
347
)
348

349
# Simulate the model
350
sim = Simulator(model, package='scipy')
351
tsim, profiles = sim.simulate(numpoints=100, integrator='odeint')
352

353
# Use simulation results to initialize discretized model
354
discretizer = TransformationFactory('dae.finite_difference') 
355
discretizer.apply_to(model, nfe=20, wrt=model.time)
356

357
# Initialize with simulation profiles
358
for t in model.time:
359
    if t in profiles:
360
        model.x[t].set_value(profiles[t]['x'])
361

362
# Now solve optimization problem
363
solver = SolverFactory('ipopt')
364
results = solver.solve(model)
365
```
366

367
### Integral Objectives and Constraints
368

369
```python
370
from pyomo.environ import *
371
from pyomo.dae import ContinuousSet, DerivativeVar, Integral
372

373
model = ConcreteModel()
374
model.time = ContinuousSet(bounds=(0, 1))
375

376
# Variables
377
model.x = Var(model.time, bounds=(0, None))
378
model.u = Var(model.time, bounds=(-2, 2))
379
model.dxdt = DerivativeVar(model.x, wrt=model.time)
380

381
# Initial condition
382
model.x[0].fix(0)
383

384
# Dynamics
385
model.ode = Constraint(
386
    model.time,
387
    rule=lambda m, t: m.dxdt[t] == m.u[t]
388
)
389

390
# Integral constraint: average control must be zero
391
model.avg_control = Integral(
392
    model.time,
393
    wrt=model.time,
394
    rule=lambda m, t: m.u[t]
395
)
396
model.avg_constraint = Constraint(expr=model.avg_control == 0)
397

398
# Integral objective: minimize integral of x^2
399
model.performance = Integral(
400
    model.time, 
401
    wrt=model.time,
402
    rule=lambda m, t: m.x[t]**2
403
)
404
model.obj = Objective(expr=model.performance, sense=minimize)
405

406
# Apply discretization
407
discretizer = TransformationFactory('dae.collocation')
408
discretizer.apply_to(model, nfe=20, ncp=2, wrt=model.time)
409
```