dmsim
simulated response of implicit/descriptor multilinear time invariant (iMTI) system to arbitrary input signals
Contents
Syntax
simout = dmsim(sys, x0) Simulated time response of autonomous iMTI system for the time fram from 0 to 10 and initial states x0.
simout = dmsim(sys, x0, ut) Simulated time response of autonomous iMTI system for the time frame of ut and initial states x0.
simout = dmsim(sys, x0, ut, u) Simulated time response of iMTI system with arbitrary input signals u (linear interpolated) at the time steps of ut and initial states x0. For continuous-time simulation use default ode-settings from Matlab.
simout = dmsim(sys, x0, ut, u, opt) Simulated time response of iMTI system with arbitrary input signals u (linear interpolated) at the time steps of ut and initial states x0, with using the ode-settings opt (defined with opt = odeset()). Only used for continuous-time simulation.
Description
The dmsim class serves as simulation class for implicit/descriptor multilinear time invariant (iMTI) models. Only differential index-1 systems are usable. The simulation will be performed automatically when creating a dmsim object by simout = dmsim(sys, x0, ut, u).
For continuous-time simulations Matlabs implementation of the variable order method ode15i() will be used. For discrete-time simulations an automated tearing approach will be used, described in [1], where the set of equations will be split into multiple subsets and explicit solvable subsets of equations will be identified, as well as implicit solvable subsets. For reordering the sparsity pattern, Matlabs implementation of the dulmage-mendelsohn decomposition dmperm() is used. The implicit solvable subsets will be solved using the Levenberg-Marquardt algorithm implementation of Matlabs nonlinear least-squares solver lsqnonlin().
Examples
Autonomous continuous-time simulation (including ode-settings): Sprott B Attractor
In [2] J. C. Sprott proposed chaotic attractor named Sprott A-S systems. The Sprott B attractor is described by the following differential equations
eq = ["dx1 = a*x2*x3";... "dx2 = x1 - b*x2";... "dx3 = c - x1*x2"];
The dmss model can be directly generated using the symbolic equations by
Ts = 0; sys = sym2dmss(eq, Ts, ["dx", "x", "u", "y"]);
where the timestep size is 0 for continuous time systems. Also it can be generated using the cpn-1 tensor format. Therefor first an empty hyCPN1() object needs to be generated, so its properties are accessable.
Ts = 0; H = hyCPN1(); H.F.stateDerivative = [1 0 0 0 0 0 0 0;... 0 1 0 0 0 0 0 0;... 0 0 1 0 0 0 0 0]; H.F.state = [0 0 0 0 1 0 0 1;... 0 0 0 1 0 1 0 1;... 0 0 0 1 0 0 0 0]; H.phi.equality = [-1 0 0 sym("a") 0 0 0 0;... 0 -1 0 0 1 -sym("b") 0 0;... 0 0 -1 0 0 0 sym("c") -1]; sys = dmss(H,Ts);
The equations can be displayed by
disp(sys.symbolicEquations)
a*x2*x3 - xp1 == 0 x1 - xp2 - b*x2 == 0 c - xp3 - x1*x2 == 0
where are the state derivatives,
are states (differential variables),
are algebriac variables and
are input signals.
For this example the parameters
params = [sym("a"), sym("b"), sym("c")];
are defined as
numParams = [0.4, 1.2, 1];
The symbolic parameter values must be replace by the numeric parameter values for the simulation.
sys = sys.replaceSymbolicParameters(params, numParams);
The Equations can be checked again by
disp(sys.symbolicEquations)
(2*x2*x3)/5 - xp1 == 0 x1 - (6*x2)/5 - xp2 == 0 1 - x1*x2 - xp3 == 0
Now the simulation can be done with creating a dmsim object with the dmss model as input, as well as initial states and a simulation time frame:
t = [0, 1e3]; x0 = [0.1, 0, 0]; opt = odeset('RelTol', 10.0^(-3),'AbsTol',10.0^(-7)); simout = dmsim(sys, x0, t, [], opt);
Seperation in 3 subset(s) with 3 subproblem(s), from which 3 are explicit solvable.
The simulation results can be accessed as the properties of the simulation object. Then the phase space can be portraited by Matlabs ploting functions
plot3(simout.x(:,1), simout.x(:,2), simout.x(:,3))
grid on

Continuous-time simulation with algebraic constraint: robertson problem
The robertson problem can be written as a system of DAEs, which gives
eq = ["0 = dx1 + 0.04*x1 - 10^4*x2*y1";... "0 = dx2 - 0.04*x1 + 10^4*x2*y1 + (3*10^7)*x2^2";... "0 = x1 + x2 + y1 - 1"]
eq = 3×1 string array "0 = dx1 + 0.04*x1 - 10^4*x2*y1" "0 = dx2 - 0.04*x1 + 10^4*x2*y1 + (3*10^7)*x2^2" "0 = x1 + x2 + y1 - 1"
As before the symbolic system of equations can be transfered into a dmss object by
Ts = 0; sys = sym2dmss(eq, Ts, ["dx", "x", "u", "y"]); disp(sys.symbolicEquations)
Reduced by 4 column(s) and 0 equation(s) with trivial Reduction due to duplications. 10000*x2*y1 - xp1 - x1/25 == 0 x1/25 - xp2 - 10000*x2*y1 - 30000000*x2*y2 == 0 1 - x2 - y1 - x1 == 0 y2 - x2 == 0
where are the state derivatives,
are states (differential variables),
are algebriac variables and
are input signals. One can see that an auxiliary variable
and an according had been added due to the quadratic term.
Simulation can be performed by creating a dmsim object
t = [0, 4e6]; x0 = [1, 0]; simout = dmsim(sys, x0, t); semilogx(simout.tsim, simout.x(:,1)) hold on semilogx(simout.tsim, simout.x(:,2)) semilogx(simout.tsim, simout.y(:,1)) grid on hold off
Seperation in 1 subset(s) with 4 subproblem(s), from which 4 are explicit solvable.

Continuous-time example with inputs: Simplified room model of an air conditioning system
HVAC systems can be modeled by dynamic energy balances (power/heat flow balances) and dynamic transport balances. The dynamics of the air temperature, enthalpy and humidity of an air conditioned room can be represented by the following power balance
powerBalance = ["0 = -m*dh1 + p1 + p2"];
transport balance of the water vapor in the air
transportBalance = ["0 = -m*dh2 + p3"];
and equation for the specific enthalpy
enthalpy = ["h1 = c_d * T1 + c_v * T1 * h2 + r0 * h2"];
where h1 is the room air enthalpy, h2 is the absolute room air humidity, dh1 and dh2 are their derivatives, T1 is the room temperature, p1 is the power transported by the air flow of the air conditioning system (either heating or cooling the room air), p2 is the heat loss to the environment and p3 is the water vapor mass flow transported by the air flow of the air conditioning system (either humidifying or dehumidifying the room air).
The equations can be converted into a continuous time iMTI model by
room = sym2dmss([powerBalance; transportBalance; enthalpy], 0, ["dh", "h", "p", "T"]);
The equations represented by the model can be viewed by
disp(room.symbolicEquations)
m*xp1 - u2 - u1 == 0 m*xp2 - u3 == 0 x1 - c_d*y1 - r0*x2 - c_v*x2*y1 == 0
where are the state derivatives,
are states (differential variables),
are algebriac variables and
are input signals.
A model for the heat and water vapor transfer by the air conditioning system can be generated by
airConditioning = ["y1 = c_d*u1*u2 + c_v*u1*u2*u3 + r0*u1*u3 - u1*u4";... "y2 = u1*u3 - u1*u5"]; airConditioning = sym2dmss(airConditioning, 0, ["dx", "x", "u", "y"]);
Reduced by 1 column(s) and 0 equation(s) with trivial Reduction due to duplications.
where y1 corresponds to the power transported by the air flow of the air conditioning system, y2 to the water vapor mass flow transported by the air flow of the air conditioning system, u1 to the supplied mass flow of air, u2 to the supply temperature, u3 to the supply humidity, u4 to the rooms enthalpy and u5 to the rooms absolute humidity.
Additionally the heat loss can be modelled by
loss = ["y1 = k*A*(u1-u2)"]; loss = sym2dmss(loss, 0, ["dx", "x", "u", "y"]);
where y1 is the heat loss, u1 is the ambient temperature and u2 is the rooms temperature.
The models can be connected in two ways. One way is using index matrices of input-state, input-algebraic and input-input connections. Therefore one of the models must be appended by the others. This can be done by
sys = append(room, airConditioning, loss);
Due to the appending the number of variables of the overall model as well as their index is increasing accordingly. The connection matrices must be produced carefully, where the first column corresponds to the input index and the second column to the state/algebraic signal index.
stateInputConnection = [7 1; 8 2]; algebraicInputConnection = [9 1; 1 2; 3 3; 2 4]; similarInputs = []; sys = sys.connect(stateInputConnection, algebraicInputConnection, [], similarInputs); disp(sys.symbolicEquations)
Reduced by 4 column(s) and 0 equation(s) with trivial Reduction due to duplications. m*xp1 - y4 - y2 == 0 m*xp2 - y3 == 0 x1 - c_d*y1 - r0*x2 - c_v*x2*y1 == 0 y2 + u1*x1 - c_d*u1*u2 - r0*u1*u3 - c_v*u1*u2*u3 == 0 y3 - u1*u3 + u1*x2 == 0 y4 + A*k*u4 - A*k*y1 == 0
Another possibility is by automatic connection by variable names, where the connection matrices will be produced automatically by identifying similar variable names of state and inputs, algebraic variables and inputs, as well as inputs and inputs.
For simulation the symbolic parameters must be replaced by numeric parameters.
symParameters = [sym("m"), sym("c_d"), sym("c_v"), sym("r0"), sym("A"), sym("k")]; numParameters = [9.03*11.33*3.4*1.225, 1006, 1860, 2501*10^3, 2*(9.03+11.33)*3.4 + 9.03*11.33, 3]; sys = sys.replaceSymbolicParameters(symParameters, numParameters);
The overall system can now be simulated with according input signals
t = [0:1:60]; x0 = [75000; 0.018]; u = [150/1.225, 18, 0.01, 32] .* ones(length(t),1); rng(100) u = u .* (0.9 + rand(size(u))/10);
where the initial state correspond to an initial room enthalpy of 75 kJ/kg and initial aboslute humidity of 18 g water vapor per kg dry air, and the (noisy) inputs correspond to an air flow of 500m³/h, supply temperature of 16°C, supply humidity of 10 g/kg and ambient temperature of 32°C. Simulation will be performed by
simout = dmsim(sys, x0, t, u); plot(simout.tsim, simout.y(:,1)) hold on plot(simout.tsim, simout.x(:,2)*1e3) legend(["room temperature in °C", "room abs. humidity in kg/kg"]) grid on axis([0 60 0 35]) hold off
Seperation in 1 subset(s) with 6 subproblem(s), from which 6 are explicit solvable.

For fast varying inputs it often makes sense to limit the maximum stepsize of the ode-solver, e.g. by opt = odeset('MaxStep', ts/2) where ts is the minimum step size between the inputs.
Discrete-time simulation using c2d
The previuous system can transformed into a discrete-time iMTI model by using the c2d() function implemented for the dmss object, using one step forward euler discretization.
Ts = 5; sys_d = sys.c2d(Ts);
System reduced algebraicly by 2 algebraic variable(s) and 2 equation(s).
Discrete-time simulation can be done similary to continuous-time simulation by creating a dmsim object. But here its important that the time vector for the input signal has the same stepsize as the stepsize of the model.
t = [0:Ts:60]; x0 = [75000; 0.018]; u = [150/1.225, 18, 0.01, 32] .* ones(length(t),1); rng(100) u = u .* (0.9 + rand(size(u))/10); sim_d = dmsim(sys_d, x0, t, u); stairs(sim_d.tsim, sim_d.y(:,1)) hold on stairs(sim_d.tsim, sim_d.x(:,2)*1e3) legend(["room temperature in °C", "room abs. humidity in kg/kg"]) grid on axis([0 60 0 35]) hold off

It can be easily seen that a stepsize of 5 seconds is not sufficient for the representation of the continuous system dynamics in discrete-time. The step size should be reduced:
Ts = 1; sys_d2 = sys_d.d2d(Ts); t = [0:Ts:60]; x0 = [75000; 0.018]; u = [150/1.225, 18, 0.01, 32] .* ones(length(t),1); rng(100) u = u .* (0.9 + rand(size(u))/10); sim_d2 = dmsim(sys_d2, x0, t, u); plot(simout.tsim, nan(length(simout.tsim), 1), 'k--', LineWidth=1.2) hold on stairs(sim_d.tsim, nan(length(sim_d.tsim), 1), 'b:', LineWidth=0.8) stairs(sim_d2.tsim, nan(length(sim_d2.tsim), 1), 'r-', LineWidth=1) plot(simout.tsim, simout.y(:,1), 'k--', LineWidth=1.2) plot(simout.tsim, simout.x(:,2)*1e3, 'k--', LineWidth=1.2) stairs(sim_d.tsim, sim_d.y(:,1), 'b:', LineWidth=0.8) stairs(sim_d.tsim, sim_d.x(:,2)*1e3, 'b:', LineWidth=0.8) stairs(sim_d2.tsim, sim_d2.y(:,1), 'r-', LineWidth=1.2) stairs(sim_d2.tsim, sim_d2.x(:,2)*1e3, 'r-', LineWidth=1.2) grid on legend(["continuous-time", "discrete-time: 5 sec.", "discrete-time: 1 sec."], "Location","northeast") hold off
System reduced algebraicly by 2 algebraic variable(s) and 2 equation(s).

References
[1] T. Warnecke and G. Lichtenberg, "Hybrid implicit multilinear simulation using difference algebraic equations reordering by sparsity patterns", 2024 10th International Conference on Control, Decision and Information Technologies (CoDIT), Valletta, Malta, 2024. In press.
[2] J. C. Sprott, “Some simple chaotic flows,” Physical Review E, vol. 50, no. 2, Aug. 1994. doi:10.1103/physreve.50.r647
See Also
dmss, mss, msim, c2d, d2c, d2d, normalize, append, connect, odeset, ode15i, lsqnonlin()
Author(s): Torben Warnecke