Azzera filtri
Azzera filtri

Find State Space Representation of Linearized non-linear System to find State Space Representation in Canonical form MATLAB

51 visualizzazioni (ultimi 30 giorni)
I have the following non-linear system and accompanying ODE:
I have the following ODE45 solution:
fun = @(t,X)odefun(X,K,C,M,F(t),resSize);
[t_ode,X_answer] = ode45(fun,tspan,X_0);
The input matrices are stiffness K(X), damping C mass M, and force F. resSize is the total number of masses in the system. To find the solution indirectly I introduced a change of variables to reduce the non-linear ODE to a first order linear ODE.
I would like to find the state space (SS) representation of my linearized non-linear system and corresponding ODE. More specifically, I would like to find the SS matrices A,B,C,D. I did not explicitly perform the change of variables. Rather, MATLAB did it for me. The ultimate goal of finding the SS representation (matrices A,B,C,D) is to then be able to find the canonical form of the non-linear ODE. I'm not sure if this can be done with a transfer function estimate of my ODE45 solution or some other MATLAB tool.

Risposta accettata

Sam Chak
Sam Chak il 17 Apr 2024
Given that no original MATLAB code for the nonlinear system is provided, I believe it would be more efficient for me to derive the linear state-space matrices intuitively by hand. If my derivation is correct, you should obtain the following symbolic state matrix :
syms m1 m2 c1 c2 c3 k1 k2 k3
%% symbolic state matrix
As = [zeros(2), eye(2);
-(k1+k2)/m1, k2/m1, -(c1+c2)/m1, 0;
k2/m2, -(k2+k3)/m2, 0, -(c2+c3)/m2]
As = 
%% symbolic input matrix
Bs = [sym('0'); sym('0'); sym('1'); sym('0')]
Bs = 
%% State-space model
A = matlabFunction(As);
B = matlabFunction(Bs);
C = eye(4);
D = 0;
m1 = 1;
m2 = 1;
c1 = 1;
c2 = 1;
c3 = 1;
k1 = 1;
k2 = 1;
k3 = 1;
sys = ss(A(c1,c2,c3,k1,k2,k3,m1,m2), B(), C, D)
sys = A = x1 x2 x3 x4 x1 0 0 1 0 x2 0 0 0 1 x3 -2 1 -2 0 x4 1 -2 0 -2 B = u1 x1 0 x2 0 x3 1 x4 0 C = x1 x2 x3 x4 y1 1 0 0 0 y2 0 1 0 0 y3 0 0 1 0 y4 0 0 0 1 D = u1 y1 0 y2 0 y3 0 y4 0 Continuous-time state-space model.
%% Step response
step(sys), grid on
  26 Commenti
Sam Chak
Sam Chak il 20 Apr 2024
The equilibrium point of the nonlinear system mentioned in my previous comment refers to the actual solution of when the system is unperturbed, specifically when . However, if the input signal u has a constant portion c, the equilibrium point may be found at solving .
As mentioned by @Torsten, the input signal u is a unbiased periodic time function, resulting in periodic oscillations observed in the simulation. Although the nonlinear system does not reach the equilibrium point in finite time, the equilibrium point still exists according to its definition. This can be likened to an undamped pendulum described by , where there are two equilibrium points, and . However, if the initial angle deviates from an equilibrium point, stable periodic oscillations will be observed.
Regarding the Jacobian matrix , which is dependent on the displacement vector , it can be computed at any given point . In this particular discussion, the Jacobian matrix mentioned is evaluated specifically at the equilibrium point. It's important to note that the forcing function is independent of the state variable x and can be treated similarly to the derivative of a constant function, .
syms x1 x2 x3 x4 m1 m2 c1 c2 c3 k1 k2 k3 gamma1 gamma2 gamma3 F
Eq1 = x3;
Eq2 = x4;
Eq3 = (F-(c1+c2)*x3-(((k1*gamma1+k2*gamma2)*x1^2+k1+k2)*x1+(-k2*gamma2*x1^2-k2)*x2))/m1;
Eq4 = ( -(c2+c3)*x4-((-k2*gamma2*x2^2-k2)*x1+((k2*gamma2+k3*gamma3)*x2^2+k2+k3)*x2))/m2;
% Jacobian matrix for a two mass network
J = jacobian([Eq1; Eq2; Eq3; Eq4], [x1, x2, x3, x4]) % Symbolic Jacobian matrix
J = 
Jonathan Frutschy
Jonathan Frutschy il 20 Apr 2024
Modificato: Jonathan Frutschy il 20 Apr 2024
%% Jacobian Matrix
J = [zeros(N), eye(N);
- M\K, - M\C]
Ok, so this is the Jacobian matrix of the system for any mass number N evaluated at the system's equilibrium point X = [x1; x_dot1; x2; x_dot2; x_3; x_dot3; x_4; x_dot4; ... x_N; x_dotN]? @Sam Chak

Accedi per commentare.

Più risposte (1)

Sam Chak
Sam Chak il 19 Apr 2024
I have conducted a simulation of the motion for 600 masses using the ode45 solver. Please refer to the comments in the code for further details, as I was unable to utilize your karray() function handle.
N = 600;
m = 1*ones(N, 1);
c = 1*ones(N+1,1);
k = 2; % <-- this one must be a scalar, or else K will be incorrect
gamma = 1;
fun = @(t, X) odefun(t, X, N, marray(m), carray(c, N), make_diagonal(X, k, gamma, N));
tspan = [0 60*pi];
X0 = zeros(2*N, 1);
[t, X] = ode45(fun, tspan, X0);
plot(t/pi, X), grid on, xlabel('t/\pi'), ylabel('\bf{x}\rm(t)'), title('System responses')
%% Multi-body Mass-Damper-Spring System
function dX = odefun(t, X, N, M, C, K)
%% Definitions
x = X(1:N); % position state vector
dx = X(N+1:2*N); % velocity state vector
%% Force vector
alpha = 0.2; % angular frequency of sin
beta = 0.15; % angular frequency of cos
famp = 0.007; % amplitude of sinusoidal wave
f = famp*sin(alpha*t)*cos(beta*t);
F = [f; zeros(175,1); f; zeros(100,1); f; zeros(74,1); f; zeros(175,1); f; zeros(71,1)];
%% Equations of Motion
ddx = M\(F - K*x - C*dx);
%% State-space model
dX = [dx; ...
ddx];
end
%% Mass matrix, M
function M = marray(m)
M = diag(m);
end
%% Damping matrix, C
function C = carray(c,N)
C = zeros(N);
for i = 1:N
C(i,i) = c(i) + c(i+1);
end
end
%% function handle for Stiffness matrix (but unused)
function K = karray(k, gamma, N)
K = @(x) make_diagonal(x, k, gamma, N);
end
%% The true Stiffness matrix, K (but k must be a scalar, which is a limitation)
function out = make_diagonal(x, k, gamma, N)
x = x(:);
x = [x(1:N); 0];
ck = circshift(k, -1);
cg = circshift(gamma, -1);
cx = circshift(x, -1);
ccx = circshift(x, 1);
d1 = -3.*ck.*cg.*cx.^2 - ck;
d2 = (k.*gamma + ck.*cg).*x.^2 + k + ck;
d3 = -3.*k.*ccx.^ 2 - k;
out = full(spdiags([d1 d2 d3], -1:1, N, N));
end
  2 Commenti
Jonathan Frutschy
Jonathan Frutschy il 20 Apr 2024
@Sam Chak This is great. I'm not sure if I was clear on how karray() works. make_diagonal() is never intended to be called. You only need to call karray() once. The only reason that my system is nonlinear is because the K matrix varies as a function of the mass displacements and velocities. The output of karray() is a function handle of @(X), that way it can be passed to odefun(), which is a function handle of @(X,t). Lets go back to the four mass example for finding the Jacobian. See below:
N = 4;
k = 2;
gamma = 1;
K_X = karray(k,gamma,N); % Function handle @(X) version of K for use in odefun()
K = K_X([zeros(N,1);zeros(N,1)]) % Find K matrix at x1(t), x2(t) = 0, x1_dot(t) = 0, and x2_dot(t) = 0
K = 4x4
4 -2 0 0 -2 4 -2 0 0 -2 4 -2 0 0 -2 4
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
%% function handle for Stiffness matrix (but unused)
function K = karray(k, gamma, N)
K = @(x) make_diagonal(x, k, gamma, N);
end
%% The true Stiffness matrix, K (but k must be a scalar, which is a limitation)
function out = make_diagonal(x, k, gamma, N)
x = x(:);
x = [x(1:N); 0];
ck = circshift(k, -1);
cg = circshift(gamma, -1);
cx = circshift(x, -1);
ccx = circshift(x, 1);
d1 = -3.*ck.*cg.*cx.^2 - ck;
d2 = (k.*gamma + ck.*cg).*x.^2 + k + ck;
d3 = -3.*k.*ccx.^ 2 - k;
out = full(spdiags([d1 d2 d3], -1:1, N, N));
end
As you can see the output for the K matrix is the same as what you got. No need to call make_diagonal(), just call karray() to get the K matrix. Does that make sense, or do you still think that my karray() function is wrong?
Jonathan Frutschy
Jonathan Frutschy il 20 Apr 2024
%% Jacobian Matrix
J = [zeros(N), eye(N);
- M\K, - M\C]
@Sam Chak How is the forcing function F(t) accounted for in this Jacobian matrix?

Accedi per commentare.

Categorie

Scopri di più su Particle & Nuclear Physics in Help Center e File Exchange

Prodotti


Release

R2024a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by