# 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