# How to simulate the response of a linear parametric varying system in Mathematica?

Consider the following system \begin{align} \dot{x}(t)&=\sum_{i=1}^{2}\rho_{i}(x(t))\left[A_{i}x(t)+B_{i}u(t)\right]\\ y(t)&=Cx(t) \end{align}

with:

a1 = {{-3, 2}, {-0.25, 1}}
a2 = {{-1.9, -0.4}, {-2.24, -4.7}}

b1 = {{0.25}, {1}}
b2 = {{-2.5}, {1}}

c = {{1, 0.5}, {0, 1}}


where

$$\rho_1(x_1(t))=\frac{1-\tanh(x_1(t))}{2},\quad\rho_2(x_1(t))=1-\rho_1(x_1(t))$$

How do I build the state-space LPV system? The response, for example, to a sinusoidal input?

I want to translate this code from MATLAB to Mathematica

A(:,:,1)=[-3 2; -0.25 1];
A(:,:,2)=[-1.9 -0.4 ; -2.24 -4.7];
B(:,:,1)=[0.25;1] ; B(:,:,2)=[-0.25; 1];
C=[1 0.05; 0 1];%
xx=[0.1; 0]; %initial conditions
mus(1,:)=[0.5 0.5 1]; %initial rhos

%%simulation loop

for k=1:tmax/Te

%% System
t(k+1)=t(k)+Te; %%time vector;

mus(1)=(1-tanh(xx(1,k)))/10  ;  % weigting functions
mus(2)=1-mus(1);

u(k+1)=10*sin(t(k)); %%input
Aa=mus(1)*A(:,:,1)+mus(2)*A(:,:,2);
Ba=mus(1)*B(:,:,1)+mus(2)*B(:,:,2);

%%Euler to solve the ODES
xx(:,k+1)=xx(:,k)+Te*(Aa*xx(:,k)+Ba*u(k) );
y(:,k+1)=C*xx(:,k+1);  %output
end

plot(t,mus,t,mus(:,1)+mus(:,2) );
plot(t,y);

• What points (in $x$ and $u$) do you wish to linearize about? – bill s Jun 14 '13 at 9:02
• What's $x_1$? And why is the argument of $\rho_i$ given by $x(t)$ in the first equation and by $t$ in the last equations? – sebhofer Jun 14 '13 at 9:02
• Do you mean a step function by $u(t)$? If you do, I think it's best to be explicit! – Vincent Tjeng Jun 14 '13 at 10:51
• $x(t)=[x_1(t)\, x_2(t)]^T$ are the states, $u$ is the input, $y(t)$ the output, and $\rho_(i)$ are the gain scheduling functions depending on the state $x_1(t)$. The system is a linear system, in fact is linear state-space parameter varying system. The idea is to simulate the response to an input $u(t)$ (for example a sinusoidal or step) similar to an standard LTI (linear time invariant system). I Matlab I made this by solving the LTI system in a for loop. But in Mathematica I don't known how to do it?. $A,B,C$ are constant matrices, $\rho_i$ changes for each step time. – user70012 Jun 14 '13 at 13:09
• This is the code for matlab – user70012 Jun 14 '13 at 13:14

StateSpaceModel will linearize the equations. The system is nonlinear. So let's look at the nonlinear solution first.

The parameters:

a1 = {{-3, 2}, {-0.25, 1}};
a2 = {{-1.9, -0.4}, {-2.24, -4.7}};
b1 = {{0.25}, {1}};
b2 = {{-2.5}, {1}};
c = {{1, 0.5}, {0, 1}};
Subscript[ρ, 1] = (1 - Tanh[Subscript[x, 1][t]])/2;
Subscript[ρ, 2] = 1 - Subscript[ρ, 1];


Set up the nonlinear equations and obtain the solution.

xx = {Subscript[x, 1][t], Subscript[x, 2][t]};
eqns = Thread[{Subscript[x, 1]'[t], Subscript[x, 2]'[t]} ==
Subscript[ρ, 1] (a1.xx + b1.{u[t]}) + Subscript[ρ, 2] (a2.xx +b2.{u[t]})];
ics = {Subscript[x, 1][0] == 0.1, Subscript[x, 2][0] == 0};
sols = NDSolve[
Join[eqns /. u[t] -> 10 Sin[t], ics], {Subscript[x, 1][t],
Subscript[x, 2][t]}, {t, 0, 10}];
p = Plot[Evaluate[c.xx /. sols], {t, 0, 10}]


The solution of the linearized system.

StateSpaceModel[eqns, xx, u[t], c.xx, t];
OutputResponse[{%, {0.1, 0}}, 10 Sin[t], {t, 0, 10}];


Compare it with the nonlinear one.

Show[p, Plot[%, {t, 0, 10}, PlotStyle -> Dashed]]


• Thanks a lot for your help. I am trying to understand all the code. Nevertheless the solution from Matlab and Mathematica are similar but quite different. Probably because in Matlab I am using Euler to solve the ODE. $\rho_i$ is a convex function, so in theory the LPV system is as the name says linear. It is possible to plot from this code the values of $\rho_i$ and also $x(t)$. Thanks for your help I am trying to change from matlab to mathematica. – user70012 Jun 14 '13 at 14:57
• First, you might try ode23 or ode45 in Matlab -- Euler approximations often go wrong. Second, this is not a linear system since $\tanh(x_1)$ is not a rational function. "Convex" and "linear" are not at all the same thing, so you should expect to see different behaviors. – bill s Jun 14 '13 at 15:00
• yes I understand, the idea of this kind of system is to have a set of LTI system [$A_ix(t)+B_iu(t)$] multiplied by a convex gain scheduling function (GSF) $\rho_i$, in this case a nonlinear function. In conclusion a set of LTIs multiplied by a nonlinear function. Where the GSF depend in the measurable (or unmeasurable) parameter varying $x(t)$. The advantage of this kind of system is that the overall representation is closer to the nonlinear behavior. Again thanks a lot for your help, I learn new things of Mathematica reading your code. – user70012 Jun 14 '13 at 15:15
• For the nonlinear system, use Plot[Evaluate[vars /. sols], {t, 0, 10}] with vars = xx for the states and vars = {Subscript[[Rho], 1], Subscript[[Rho], 2]} for the parameters. For StateSpaceModel use StateResponse for the states, and change the outputs to obtain the parameters. Regarding nonlinear, I used it in the sense of not being linear (superposition and homogeneity fails.) – Suba Thomas Jun 14 '13 at 15:21
• Your right. People in control also call this kind of systems as Takagi-sugeno. Thanks for your help. – user70012 Jun 14 '13 at 15:27

## Update

As of Version 10 nonlinear control systems can be addressed by NonlinearStateSpaceModel. With this the answer given by @Suba Thomas can be obtained by:

a1 = {{-3, 2}, {-0.25, 1}};
a2 = {{-1.9, -0.4}, {-2.24, -4.7}};
b1 = {{0.25}, {1}};
b2 = {{-2.5}, {1}};
c = {{1, 0.5}, {0, 1}};
ρ[1] = (1 - Tanh[x[1][t]])/2;
ρ[2] = 1 - ρ[1];

xx = { x[1][t], x[2][t] };
eqns = Thread[{x[1]'[t], x[2]'[t]} == ρ[1] (a1.xx + b1.{u[t]})
+ ρ[2] (a2.xx + b2.{u[t]})];
ics = {0.1, 0}; (* initial conditions for x[1][0], x[2][0] *)

(* Nonlinear State Space Model *)
nssm = NonlinearStateSpaceModel[ eqns, xx, u[t], c.xx, t ];
nssmResp = OutputResponse[ { nssm, ics }}, 10 Sin[t], {t, 0, 10} ];
nssmStyle = Sequence[ Directive[Blue], Directive[Red] ];

(* Linearized State Space Model *)
ssm = StateSpaceModel[eqns, xx, u[t], c.xx, t];
ssmResp = OutputResponse[ { ssm, ics }, 10 Sin[t], {t, 0, 10} ];
ssmStyle = Sequence[Directive[ Dashed, Blue], Directive[Dashed, Red]];

(* Plotting the Responses *)
legend = LineLegend[
{ Blue, Red, Directive[Dashed, Blue], Directive[Dashed, Red] },
{"nonlinear g(x1)", "nonlinear g(x2)", "linearized g(x1)", "linearized g(x2)"}
];
Plot[ {ssmResp, nssmResp}, {t, 0, 10},
PlotStyle -> {ssmStyle, nssmStyle},
PlotLegends -> legend
]


• Note that I shun Subscript as the devil will holy water (repeating myself)... – gwr May 25 '16 at 15:14
• Yes, this is now the more natural way to obtain the solution after NonlinearStateSpaceModel came onto the scene. – Suba Thomas May 26 '16 at 14:51
• @SubaThomas I must admit that (not coming from engineering) I am still struggling with the question when to use something like ParametricNDSolve or the block based modeling of NonlinearStateSpaceModel. There seem to be ups and downs for both: WhenEvent, speed of simulation when varying a lot of parameters, state resetting for maybe Bayesian estimation ... Maybe that is a question worth asking in itself? – gwr May 26 '16 at 15:03