Inconsistent result: Integration of step response vs direct evaluation

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]

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.

Answer 1

[Sam Chak]

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]

Thanks, but do you have an idea how the integral would actually work ? Here is the snippet to equation 84 below, the plot is blank. If no response, i may have to create a new topic. I am interested in the integral beacuse i am looking at other response function as well. So, i want to know how the integral with tau works. thanks:

clear all; close all; clc

syms tau t real

A =  [-1 1; -4 -4];
B = [0; 4];
K = 10;


y(t) = int(exp(A*(t-tau))*B*K,tau,0,t,'Hold',true);
fplot(y(t), [0 2], 'c');

[Torsten]

clear all; close all; clc

syms tau t real

A =  [-1 1; -4 -4];
B = [0; 4];
K = 10;


y(t) = int(expm(A*(t-tau))*B*K,tau,0,t);
fplot(y(t), [0 2], 'c');

[Sam Chak]

@Torsten, Thank you!

[Joe]

Thanks Torsten and Sam Chak. But do you have an idea how to get the actual values from the inetgration so they can be saved in a .MAT file for instance because I see the worskpace has symbolic variables not actual numbers. Thank you

[Torsten]

syms tau t real
A =  [-1 1; -4 -4];
B = [0; 4];
K = 10;
y = int(expm(A*(t-tau))*B*K,tau,0,t);
T = 0:0.1:2;
Y = double(subs(y,t,T))

Y = 2×21

0 0.1694 0.5742 1.0965 1.6574 2.2067 2.7145 3.1654 3.5537 3.8796 4.1474 4.3631 4.5338 4.6664 4.7678 4.8438 4.8996 4.9398 4.9679 4.9870 4.9993 0 3.2755 5.3700 6.6173 7.2726 7.5280 7.5248 7.3649 7.1201 6.8390 6.5535 6.2832 6.0390 5.8259 5.6450 5.4949 5.3730 5.2757 5.1994 5.1407 5.0963

hold on
plot(T,Y(1,:),'r')
plot(T,Y(2,:),'b')
hold off
grid on

[Joe]

Thank yoo Torsten

[Joe]

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)
Xstr = mupadmex('symobj::double', S.s, 0);

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]

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

Answer 2

[Sam Chak]

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