Add consistent initialization for unit operation models (or residuals)

[daklauss]

In GitLab by @schmoelder on Jul 4, 2024, 17:30

To run the simulation, we also need a method to consistently initialize the system for given initial values.

Here is an example for how to initialize the CSTR:

# differential equations
d (c*V) / dt = Q_in * c_in - Q_out * c
d V / dt = Q_in - Q_out

# residual (c)
0 = c_dot * V0 + V_dot * c - Q_in * c_in + Q_out * c

# residual (V)
0 = V_dot - Q_in + Q_out

# Problem: For given initial values (c0, V0), what are the initial values for V_dot0, c_dot0? => Solve residual equations.
# Here, we need to start with V0 since it's also required for c_dot0
V_dot0 = Q_in - Q_out
c_dot0 = - V_dot0/V0 * c0 + Q_in/V0 * c_in - Q_out/V0 * c0

A method to compute the initial values needs to be added for all residual implementations.

Open question: Should we maybe add a class for the residual?

[daklauss]

Consistent initialization with initial values is dependend on the Unit Operation

• Spliting UnitOperationBase.initialize() in initialize_state(), initialize_state_derivative
• Calling unit_operation.initialize() initializes the state
• User then provides y,
• ydot is calculated by calling initialize_state_derivative
• Known iv depends on the Unit Operation and coupling

[daklauss]

Create a initialize_state_derivative for the system. This function uses scipy.optimize.fsolve to calculate the remaining unknown derivatives that are dependend on the connection, that can't be calculated directly. Until fsolve is called, it iterates over each unit operation to call their own initialize_state_derivative that tries to calculate their initial values and gives the fsolve an initial value.

[daklauss]

It is needed to add to calculate the initial values for each unit operation:
Inlet should be as follows if you consider cadet:

# residual c:
t_poly = (1, t, t², t³)
0 = c_poly @ t_poly - c

# so c_dot should be
t_poly = (0, 1, 2t, 3t²)
c_dot = c_poly @ t_poly

[schmoelder]

@daklauss Just out of curiosity, does this now work for all cases or are there still some situations, where we need to add some functionality (e.g. cyclic systems)?

[daklauss]

It does not work for all cases. The thing is, this is not an issue of specific operation models but rather an issue of the whole system. I would rather open an own issue for that and rename it accordingly.