# Inconsistent result: Integration of step response vs direct evaluation

1 view (last 30 days)
Joe on 19 Aug 2023
Edited: Torsten on 22 Aug 2023
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.
Now, the example problem:::
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)
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);
Error using vertcat
Dimensions of arrays being concatenated are not consistent.

Error in solution>testode1 (line 35)
dy = [exp(-ii*t); y(1)];

Error in odearguments (line 92)
f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 107)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);
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
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.

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
Joe on 21 Aug 2023
Edited: Joe on 21 Aug 2023
I just stunbled on a problem, sorry man, Suppose I want to mutiply B by a function of t. This funcction has an if /elseif and else for loop defined below. The system falls apart because t is not defined.
As an additional question, what if the value of y(t) vector is needed for evaluation in the Function during the integration process, how would this go..... For example instead of (if 5+double(subs(t))>= 5.5) in the funcC function, replace it with if (y(2)>= 5.5). I am curious about this. Some evaluation may need some values in a vector for evaluating a piecewise function for the next step. Thanks
syms tau t real
A = [-1 1; -4 -4];
B(t) = [0; 4]*funcC(t);
Error using symengine
Unable to convert expression containing symbolic variables into double array. Apply 'subs' function first to substitute values for variables.

Error in sym/double (line 729)

Error in solution>funcC (line 18)
if 5+double(subs(t))>= 5.5
K = 10;
y = int(expm(A*(t-tau))*B(t)*K,tau,0,t);
T = 0:0.1:2;
Y = double(subs(y,t,T))
hold on
plot(T,Y(1,:),'r')
plot(T,Y(2,:),'b')
hold off
grid on
function [dx] = funcC(t)
if 5+double(subs(t))>= 5.5
dx = 2;
elseif 5+double(subs(t))<= 3.5
dx =35;
else
dx = 0.55;
end
end
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)
B =
K = 10;
y = int(expm(A*(t-tau))*B*K,tau,0,t);
T = 0:0.01:2;
Y = double(subs(y,t,T))
Y = 2×201
0 0.0011 0.0043 0.0094 0.0165 0.0253 0.0358 0.0480 0.0616 0.0767 0.0931 0.1109 0.1298 0.1498 0.1709 0.1929 0.2159 0.2397 0.2644 0.2897 0.3158 0.3425 0.3697 0.3975 0.4258 0.4545 0.4836 0.5130 0.5428 0.5728 0 0.2156 0.4227 0.6216 0.8124 0.9953 1.1708 1.3389 1.4999 1.6540 1.8015 1.9425 2.0773 2.2060 2.3289 2.4461 2.5578 2.6643 2.7656 2.8619 2.9535 3.0404 3.1229 3.2011 3.2751 3.3451 3.4112 3.4735 3.5323 3.5876
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

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:
.
Suppose the system input u(t) is given by
,
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
The solution for the reference input is provided by the following formula:
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
Phi =
% Solution of Nonhomogeneous Linear State Equations
x(t) = Phi*x0 + int(expm(Ac*(t - tau))*B*N*r, tau, 0, t)
x(t) =
y(t) = C*x(t)
y(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

### Categories

Find more on Symbolic Math Toolbox in Help Center and File Exchange

R2022b

### Community Treasure Hunt

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

Start Hunting!