2
$\begingroup$

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); 
$\endgroup$
7
  • $\begingroup$ What points (in $x$ and $u$) do you wish to linearize about? $\endgroup$
    – bill s
    Jun 14, 2013 at 9:02
  • $\begingroup$ 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? $\endgroup$
    – sebhofer
    Jun 14, 2013 at 9:02
  • $\begingroup$ Do you mean a step function by $u(t)$? If you do, I think it's best to be explicit! $\endgroup$ Jun 14, 2013 at 10:51
  • $\begingroup$ $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. $\endgroup$
    – user70012
    Jun 14, 2013 at 13:09
  • $\begingroup$ This is the code for matlab $\endgroup$
    – user70012
    Jun 14, 2013 at 13:14

2 Answers 2

7
$\begingroup$

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}]

enter image description here

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]]

enter image description here

$\endgroup$
5
  • $\begingroup$ 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. $\endgroup$
    – user70012
    Jun 14, 2013 at 14:57
  • $\begingroup$ 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. $\endgroup$
    – bill s
    Jun 14, 2013 at 15:00
  • $\begingroup$ 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. $\endgroup$
    – user70012
    Jun 14, 2013 at 15:15
  • $\begingroup$ 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.) $\endgroup$ Jun 14, 2013 at 15:21
  • $\begingroup$ Your right. People in control also call this kind of systems as Takagi-sugeno. Thanks for your help. $\endgroup$
    – user70012
    Jun 14, 2013 at 15:27
2
$\begingroup$

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
]

OutputResponses

$\endgroup$
3
  • $\begingroup$ Note that I shun Subscript as the devil will holy water (repeating myself)... $\endgroup$
    – gwr
    May 25, 2016 at 15:14
  • $\begingroup$ Yes, this is now the more natural way to obtain the solution after NonlinearStateSpaceModel came onto the scene. $\endgroup$ May 26, 2016 at 14:51
  • $\begingroup$ @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? $\endgroup$
    – gwr
    May 26, 2016 at 15:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.