Inconsistent result: Integration of step response vs direct evaluation
1 view (last 30 days)
Show older comments
I am following a work from MIT on state space -: step response. I will present 2 pages of the text and then the example problem.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1460037/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1460042/image.png)
Now, the example problem:::
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1460047/image.png)
My question and problem next:
I have done equation (vi) using a for loop below and the result is correct as the text, see below:
tt=0:1e-3:2;
for q=1:length(tt)
t=tt(q);
pp(q,:)= 5.0 + exp(-2.5*t)*(-5*cos(1.323*t) + 20.8*sin(1.323*t));
end
plot(tt,pp)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1460057/image.png)
However, when I tried to integrate the actual equation in eqn(84), problem happens. I think it is because the A matrix in eqn 84 has no eigen values like eqn vi. So, I have used the eigen matrix (2x2) used in eqn(v) in the last image as a replacement. However, it is not right. Anyone knows how the actual interation would work can help,,,, thank you My code below. It shows an error due to a 2x2 from exp(-ii*t) in the ode function.Thank you very much in advance
clear all; close all; clc
global ii A B K
% % % cam = [-0.3125-0.1654*j -0.3125+0.1654*j];
% % % mbk = [-20 + 22.678*j; -20 - 22.678*j];
% % % cabk = [-5.0];
ii = [-2.5+1.3223*j 0; 0 -2.5-1.3223*j];
A = [-1 1; -4 -4];
B = [0; 4];
K = 10;
dt = .0002; %Defining Time Step:
t = 0:dt:2;% Defining time vector.
y0 = [0; 0]; %initial vel and disp [vel disp] %velocity.
[tsol, yvectorl] = ode45(@testode1,t,y0);
yvectorl = exp(-ii*0)*yvectorl*B*K;
%%plot
figure()
plot(tsol,yvectorl(:,2)), title("Displacement") %Disp.(2)
figure()
plot(tsol,yvectorl(:,1)), title("Velocity") %Vel. (1
function dy = testode1(t,y)
global ii A B K
dy = [exp(-ii*t); y(1)]; %%%error due to a 2x2 from exp(-ii*t) here
end
1 Comment
Paul
on 19 Aug 2023
I realize that I'm only seeing a snippet of the material ....
u_s(t) shouldn't be in equation (86) because it's already been assumed that u_s(t>=0) = 1 in the derivation of equation (86).
The statement immediately preceding equation (87) is incorrect. The system will only have a steady state response if a) all of the eigenvalues of A are in the open left half plane, or b) any eigenvalues of A in the closed right half plane are cancelled by zeros, as could be the case if A,B,C,D is a non-minimal realization.
Accepted Answer
Sam Chak
on 19 Aug 2023
Hi @Joe
The ode45() function can be implemented this way:
tspan = 0:0.01:2;
x0 = [0 0]; % initial condition
[t, x] = ode45(@odefcn, tspan, x0);
plot(t, x(:, 2), 'linewidth', 1.5); grid on;
hold on
y = 5.0 + exp(-2.5*t).*(-5*cos(1.323*t) + 20.8*sin(1.323*t));
plot(t, y)
hold off
legend('numerical sol', 'analytical sol')
xlabel('t'), ylabel('x(t)')
function dxdt = odefcn(t, x)
A = [-1 1; -4 -4];
B = [0; 4];
C = [0 1];
D = 0;
K = 10;
u = K;
dxdt = A*x + B*u;
end
8 Comments
Torsten
on 21 Aug 2023
Edited: Torsten
on 22 Aug 2023
Use "piecewise" to define B - this is a symbolic function.
You should compare with the numerical solution using ODE45 because of the jump discontinuity.
syms tau t real
A = [-1 1; -4 -4];
B = [0; 4]*piecewise(tau<=-1.5,35,(tau>-1.5)&(tau<0.5),0.55,tau>=0.5,2)
K = 10;
y = int(expm(A*(t-tau))*B*K,tau,0,t);
T = 0:0.01:2;
Y = double(subs(y,t,T))
hold on
plot(T,Y(1,:),'r')
plot(T,Y(2,:),'b')
x0 = [0 0]; % initial condition
tspan = 0:0.01:0.5;
iflag = 1;
[t1, x1] = ode45(@(t,x)odefcn(t,x,iflag), tspan, x0);
x0 = x1(end,:);
tspan = 0.5:0.01:2;
iflag = 2;
[t2, x2] = ode45(@(t,x)odefcn(t,x,iflag), tspan, x0);
t = [t1;t2(2:end)];
x = [x1;x2(2:end,:)];
plot(t, x(:,1),'r')
plot(t, x(:,2),'b')
hold off
grid on
function dxdt = odefcn(t, x, iflag)
A = [-1 1; -4 -4];
if iflag == 1
B = [0; 4]*0.55;
elseif iflag == 2
B = [0; 4]*2;
end
C = [0 1];
D = 0;
K = 10;
u = K;
dxdt = A*x + B*u;
end
More Answers (1)
Sam Chak
on 19 Aug 2023
Hi @Joe
For completeness, I am including an example of the full-state variable feedback simulation using the symbolic matrix exponential (expm) approach. Consider the state-space system:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1460232/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1460237/image.png)
Suppose the system input u(t) is given by
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1460242/image.png)
where r(t) is the reference input, and N is an input scaling factor designed to eliminate the steady-state tracking error for a step input. Consequently, the closed-loop system equation takes the form
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1460247/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1460252/image.png)
The solution for the reference input is provided by the following formula:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1460257/image.png)
syms t tau real
A = [0 1;
0 0]; % state matrix
B = [0;
1]; % input matrix
C = [1 0]; % output matrix
K = [1 2]; % feedback gain matrix
Ac = A - B*K; % closed-loop state matrix
r = 1; % reference input, unit step input, in this example
N = 1; % input scaling factor
x0 = [0; 0]; % initial condition
Phi = expm(Ac*t) % state transition matrix
% Solution of Nonhomogeneous Linear State Equations
x(t) = Phi*x0 + int(expm(Ac*(t - tau))*B*N*r, tau, 0, t)
y(t) = C*x(t)
figure(1)
fplot(y(t), [0 10]), grid on
title('Unit Step Response')
xlabel({'$t$'}, 'interpreter', 'latex', 'fontsize', 16)
ylabel({'$y(t)$'}, 'interpreter', 'latex', 'fontsize', 16)
Furthermore, it's worth noting that the unit step response displayed in Figure 1 yields the same response as shown in Figure 2 below, which is generated using the step() function.
figure(2)
sys = ss(Ac, B*N, C, 0);
step(sys, 10), grid on
0 Comments
See Also
Categories
Find more on Symbolic Math Toolbox in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!