Model Differential Algebraic Equations

Overview of Robertson Reaction Example

Robertson[1]created a system of autocatalytic chemical reactions to test and compare numerical solvers for stiff systems. The reactions, rate constants (k), and reaction rates (V) for the system are given as follows:

$\begin{matrix} A & \overset{k_{1}}{\to} & B & k_{1} = 0.04 & V_{1} = k_{1} [A] \\ B + B & \overset{k_{2}}{\to} & C + B & k_{2} = 3 \cdot 10^{7} & V_{2} = k_{2} [B] [B] \\ B + C & \overset{k_{3}}{\to} & A + C & k_{3} = 1 \cdot 10^{4} & V_{3} = k_{3} [B] [C] \end{matrix}$

Because there are large differences between the reaction rates, the numerical solvers see the differential equations as stiff. For stiff differential equations, some numerical solvers cannot converge on a solution unless the step size is extremely small. If the step size is extremely small, the simulation time can be unacceptably long. In this case, you need to use a numerical solver designed to solve stiff equations.

Simulink Model from ODE Equations

A system of ordinary differential equations (ODE) has the following characteristics:

All of the equations are ordinary differential equations.
Each equation is the derivative of a dependent variable with respect to one independent variable, usually time.
The number of equations is equal to the number of dependent variables in the system.

Using the reaction rates, you can create a set of differential equations describing the rate of change for each chemical species. Since there are three species, there are three differential equations in the mathematical model.

$\begin{array}{l} A^{'} = - 0.04 A + 1 \cdot 10^{4} B C \\ B^{'} = 0.04 A - 1 \cdot 10^{4} B C - 3 \cdot 10^{7} B^{2} \\ C^{'} = 3 \cdot 10^{7} B^{2} \end{array}$

Initial conditions: $A = 1$ , $B = 0$ , and $C = 0$ .

Build the Model

Create a model, or open the modelex_hb1ode.

Add three Integrator blocks to your model. Label the inputsA',B', andC', and the outputsA,B, andCrespectively.
Add Sum, Product, and Gain blocks to solve each differential variable. For example, to model the signalC',
1. Add a Math Function block and connect the input to signalB. Set theFunctionparameter tosquare.
2. Connect the output from the Math Function block to a Gain block. Set theGainparameter to3e7.
3. 继续添加剩余的微分equation terms to your model.
Model the initial condition ofAby setting theInitial conditionparameter for theAIntegrator block to1.
Add Out blocks to save the signalsA,B, andCto the MATLAB variableyout.

Simulate the Model

Create a script that uses thesimcommand to simulate your model. This script saves the simulation results in the MATLAB variableyout. Since the simulation has a long time interval andBinitially changes very fast, plotting values on a logarithmic scale helps to visually compare the results. Also, since the value ofBis small relative to the values ofAandC, multiplyBby $1 \cdot 10^{4}$ before plotting the values.

Enter the following statements in a MATLAB^®script. If you built your own model, replaceex_hblodewith the name of your model.

sim('ex_hb1ode') yout(:,2) = 1e4*yout(:,2); figure; semilogx(tout,yout); xlabel('Time'); ylabel('Concentration'); title('Robertson Reactions Modeled with ODEs')

From the Simulink^®Editor, on theModelingtab, clickModel Settings.
— In the Solver pane, set theStop timeto4e5and theSolvertoode15s (stiff/NDF).
— In the Data Import pane, select theTimeandOutputcheck boxes.
Run the script. Observe that all ofAis converted toC.

Simulink Model from DAE Equations

A system of differential algebraic equations (DAE) has the following characteristics:

It contains both ordinary differential equations and algebraic equations. Algebraic equations do not have any derivatives.
Only some of the equations are differential equations defining the derivatives of some of the dependent variables. The other dependent variables are defined with algebraic equations.
The number of equations is equal to the number of dependent variables in the system.

Some systems contain constraints due to conservation laws, such as conservation of mass and energy. If you set the initial concentrations to $A = 1$ , $B = 0$ , and $C = 0$ , the total concentration of the three species is always equal to $1$ since $A + B + C = 1$ . You can replace the differential equation for $C^{'}$ with the following algebraic equation to create a set of differential algebraic equations (DAEs).

$C = 1 - A - B$

The differential variablesAandBuniquely determine the algebraic variableC.

$\begin{array}{l} A^{'} = - 0.04 A + 1 \cdot 10^{4} B C \\ B^{'} = 0.04 A - 1 \cdot 10^{4} B C - 3 \cdot 10^{7} B^{2} \\ C = 1 - A - B \end{array}$

Initial conditions: $A = 1$ and $B = 0$ .

Build the Model

Make these changes to your model or to the modelex_hb1ode, or open the modelex_hb1dae.

Delete the Integrator block for calculatingC.
Add a Sum block and set theList of signsparameter to +– –.
Connect the signalsAandBto the minus inputs of the Sum block.
Model the initial concentration ofAwith a Constant block connected to the plus input of the Sum block. Set theConstant valueparameter to1.
Connect the output of the Sum block to the branched line connected to the Product and Out blocks.

Simulate the Model

Create a script that uses thesimcommand to simulate your model.

Enter the following statements in a MATLAB script. If you built your own model, replaceex_hbldaewith the name of your model.

sim('ex_hb1dae') yout(:,2) = 1e4*yout(:,2); figure; semilogx(tout,yout); xlabel('Time'); ylabel('Concentration'); title('Robertson Reactions Modeled with DAEs')

From the Simulink Editor, on theModelingtab, clickModel Settings.
— In the Solver pane, set theStop timeto4e5and theSolvertoode15s (stiff/NDF).
— In the Data Import pane, select theTimeandOutputcheck boxes.
Run the script. The simulation results when you use an algebraic equation are the same as for the model simulation using only differential equations.

Simulink Model from DAE Equations Using Algebraic Constraint Block

You can replace the differential equation for $C^{'}$ with an algebraic equation modeled using an Algebraic Constraint block and a Sum block. The Algebraic Constraint block constrains its input signal F(z) to zero and outputs an algebraic state z. In other words, the block output is a value needed to produce a zero at the input. Use the following algebraic equation for input to the block.

$0 = A + B + C - 1$

The differential variablesAandBuniquely determine the algebraic variableC.

$\begin{array}{l} A^{'} = - 0.04 A + 1 \cdot 10^{4} B C \\ B^{'} = 0.04 A - 1 \cdot 10^{4} B C - 3 \cdot 10^{7} B^{2} \\ C = 1 - A - B \end{array}$

Initial conditions: $A = 1$ , $B = 0$ , and $C = 1 \cdot 10^{- 3}$ .

Build the Model

Make these changes to your model or to the modelex_hb1ode, or open the modelex_hb1dae_acb.

Delete the Integrator block for calculatingC.
Add an Algebraic Constraint block. Set theInitial guessparameter to1e-3.
Add a Sum block. Set theList of signsparameter to –+++.
Connect the signalsAandBto plus inputs of the Sum block.
Model the initial concentration ofAwith a Constant block connected to the minus input of the Sum block. Set theConstant valueparameter to1.
Connect the output of the Algebraic Constraint block to the branched line connected to the Product and Out block inputs.
Create a branch line from the output of the Algebraic Constraint block to the final plus input of the Sum block.

Simulate the Model

Create a script that uses thesimcommand to simulate your model.

Enter the following statements in a MATLAB script. If you built your own model, replaceex_hbl_acbwith the name of your model.

sim('ex_hb1dae_acb') yout(:,2) = 1e4*yout(:,2); figure; semilogx(tout,yout); xlabel('Time'); ylabel('Concentration'); title('Robertson Reactions Modeled with DAEs and Algebraic Constraint Block')

From the Simulink Editor, on theModelingtab, clickModel Settings.
— In the Solver pane, set theStop timeto4e5and theSolvertoode15s (stiff/NDF).
— In the Data Import pane, select theTimeandOutputcheck boxes.
Run the script. The simulation results when you use an Algebraic Constraint block are the same as for the model simulation using only differential equations.

Using an Algebraic Constraint block creates an algebraic loop in a model, If you set theAlgebraic Loopparameter towarning(on theModelingtab, clickModel Settings, then selectDiagnostics), the following message displays in the Diagnostic Viewer during simulation.

For this model, the algebraic loop solver was able to find a solution for the simulation, but algebraic loops don’t always have a solution, and they are not supported for code generation. For more information about algebraic loops in Simulink models and how to remove them, seeAlgebraic Loop Concepts.

References

[1] Robertson, H. H. “The solution of a set of reaction rate equations.”Numerical Analysis: An Introduction(J. Walsh ed.). London, England:Academic Press, 1966, pp. 178–182.

Related Examples

Best-Form Mathematical Models