Linearization of the ODE system: Problems

Question

I have summarized the issues covered in the topics:

https://mathematica.stackexchange.com/questions/253133/linearization-of-ode-without-an-equilibrium

https://mathematica.stackexchange.com/questions/253144/i-ask-for-help-with-commands-transferfunctionmodel-statespacemodel

https://mathematica.stackexchange.com/questions/253092/plot3d-whenevent-ndsolve

---

System of ODE:

$\begin{cases} \ddot{s}+\omega^2 \cdot s=0 \\ \dot{x}=G+u \\ \dot{z}+z=\frac{d}{dt}(-(x+s-1)^2) \\ \dot{G}+G=z \cdot s \\ \end{cases}$

Problem:

All codes in Mathematica:

Given the ODE system:

Clear["Derivative"]; ClearAll["Global`*"];

pars = {\[Alpha] = 1, \[Omega] = 2 Pi 2};

asys = AffineStateSpaceModel[{s''[t] + \[Omega]^2 s[t] == 0, x'[t] == G[t] + u[t], z'[t] + z[t] == D[-(x[t] + s[t] - 1)^2, t], G'[t] + G[t] == z[t] s[t]}, {{s[t], 0}, {s'[t], \[Alpha] \[Omega]}, {x[t], 0}, {z[t], 0}, {G[t], 0}}, {u[t]}, {x[t], s[t], z[t], G[t]}, t];

---

Calculating the required output signal:

or = OutputResponse[asys, 0, {t, 0, 25}];

Plot[{or[[1]], 1}, {t, 0, 15}, PlotRange -> Full];

If we do not take small oscillations (we do not need them) contained in the trajectory, then it is clear that this oscillatory transient process can be roughly described by a quasilinear law, which means that a certain transfer function corresponds to it.

---

Drawing on ideas from:

https://mathematica.stackexchange.com/questions/138852/linearization-of-a-nonlinear-system

As well as:

https://mathematica.stackexchange.com/questions/246111/system-of-ode-rightarrow-affine-state-space-rightarrow-system-of-ode-in-ca

I am trying to linearize the system at the equilibrium point using the Jacobi matrix. For some reason, one of the balance points is not searched for, so I set it manually:

ClearAll[stateEqs]; stateEqs[Verbatim[AffineStateSpaceModel][{a_, b_, c_, d_}, x0_, u0_, y_, t_]] := Module[{u, x}, x = Replace[x0, {xx_, x1_} :> xx, 1]; u = Replace[u0, {uu_, u1_} :> uu, 1]; D[x, t] == a + b.u];
 
newode = asys // stateEqs // Thread; newode // ExpandAll // MapAt[Simplify, #, {All, 2}] & // Column // Expand;


eqns = ReplaceAll[newode /. Thread[Array[Subscript[\[FormalX], #] &, 5] -> {ds, s, x, z, G}] // DeleteCases[True], u[t] -> 0] // Expand;

eq = Solve[{s[t], -16 \[Pi]^2 ds[t], G[t], 2 G[t] - 2 ds[t] G[t] + 2 s[t] - 2 ds[t] s[t] - 2 G[t] x[t] - 2 s[t] x[t] - z[t], -G[t] + ds[t] z[t]} == {0, 0, 0, 0, 0}, {ds[t], s[t], x[t], z[t], G[t]}, Reals] // Simplify; (*Equilibrium Points*)

eq = {{ds[t] -> 0, s[t] -> 0, x[t] -> 0, z[t] -> 0, G[t] -> 0}}(*Handly equilibrium point*)

---

Here I calculate the Jacobi matrix:

j = D[{s[t], -16 \[Pi]^2 ds[t], G[t], 2 G[t] - 2 ds[t] G[t] + 2 s[t] - 2 ds[t] s[t] - 2 G[t] x[t] - 2 s[t] x[t] - z[t], -G[t] + ds[t] z[t]}, {{ds[t], s[t], x[t], z[t], G[t]}}];(*Jacobi Matrix*)

J = N[j /. eq[[1]]] // Simplify;

eqns = Thread[{ds'[t], s'[t], x'[t], z'[t], G'[t]} == J.{ds[t], s[t], x[t], z[t], G[t]} + {0, 0, u[t], 0, 0}];(*Collect ODE's*)

asys = AffineStateSpaceModel[eqns, {{s[t], 0}, {s'[t], \[Alpha] \[Omega]}, {x[t], 0}, {z[t], 0}, {G[t], 0}}, {u[t]}, {x[t]}, t](*Linearized Affine State-Space*)

---

When calculating the output of a "linearized" system, I get sheer nonsense. Either I made a mistake in determining the stationary point, or something else.

or2 = OutputResponse[asys, 0, {t, 0, 25}];

Plot[{or2}, {t, 0, 25}, PlotRange -> Full, PlotPoints -> 200];

Resume:

Dear specialists.With all the obvious similarity of the transient processes in the original system, my attempts to describe it with a linear state-space repeatedly fail. Help me figure out this problem. Topic will be useful to others, and also give me a powerful impetus to understand the problem and master the tools. I would be grateful to everyone for help and advice .

Answer 1

I will first write your system of differential equations as a system of first order differential equations and group some terms. It is, $$\begin{aligned} \dot{s} &= q,\\ \dot{q} &= -\omega^2\, s,\\ \dot{x} &= G + u,\\ \dot{z} &= -z - 2(x + s - 1)(G - q) - 2(x + s - 1)\,u,\\ \dot{G} &= -G + z\, s \end{aligned}$$ Mathematica does produce the right model from your equations and, if you run asys through StateSpaceModel, it produces the right linearization too:

You do not need to compute it manually as you try to. However, as Kwin points out, the correct linearization has eigenvalues on the imaginary axis. This means that the behaviour of the linearized system does not mimic the local behaviour of the original nonlinear system. The eigenvalues, by the way, are $\pm 4 \pi\,i,$ $-1,$ $-1,$ and $0.$ As a result, there is no use in acquiring the output response of the linearized system since it will tell you nothing useful; any conclusion you draw from it would be unfounded.

Even if that was not the case, your output response is expected. You computed the output response of a linear state-space model with zero input. The output is zero in that case. The output response computes the output $y(t)$ given an input $u(t)$ assuming the initial condition is at the equilibrium point. Compute the step response (pass $1$ instead of $0$) and you'll see something. If you instead compute the state response (using StateResponse) with a non-zero initial condition for the linear system, then you will observe non-trivial behaviour.

All of that said, you can still perform linear control design and stabilize the system locally if you really want to. But that should address your question about the output of the linearized system.