Evaluate Settling Time,Constant Time of my System

Hello everyone,i am studying my system,it'about an LQR control of a double pendulum around 0°angle equilibrium.I have to discuss about how poles of my system moves,changing Q and R values and i have to discuss about how Q and R high values could make my system faster but potentially create oscillations and in opposite case,slower...so i wanted to analize the Settling Time of my LQR linear system states plots to show that..How could i evaluate from my plots the settling time.I saw something around stepinfo but that was the response of system from step input in my case there is no input..i tried something to see when values where higher than 63% of final state value but it didnt work..can someone please help me.thanks

m0  = 1.5;
m1  = 0.5;
m2  = 0.75;
L1  = 0.5;
L2  = 0.75;
g   = 9.8;
d1  = m0 + m1 + m2;
d2  = (0.5*m1 + m2)*L1;
d3  = 0.5*m2*L2;
d4  = (1/3*m1 + m2)*L1^2;
d5  = 0.5*m2*L1*L2;
d6  = 1/3*m2*L2^2;
f1  = (0.5*m1 + m2)*L1*g;
f2  = 0.5*m2*L2*g;
%% Symbolic matrices in the nonlinear dynamics
syms Theta0 Theta1 Theta2 DTheta0 DTheta1 DTheta2;
x_  = transpose([Theta0 Theta1 Theta2 DTheta0 DTheta1 DTheta2]);
D(x_)= [d1            d2*cos(x_(2))       d3*cos(x_(3));
        d2*cos(x_(2)) d4                  d5*cos(x_(2)-x_(3));
        d3*cos(x_(3)) d5*cos(x_(2)-x_(3)) d6];
G(x_)= [0; 
       -f1*sin(x_(2)); 
       -f2*sin(x_(3))];
C(x_)=[0 -d2*sin(x_(2))*x_(5)      -d3*sin(x_(3))*x_(6);
       0  0                         d5*sin(x_(2)-x_(3))*x_(6);
       0 -d5*sin(x_(2)-x_(3))*x_(5) 0];
H    = [1; 
        0; 
        0];
%% matrici A e B sistema linearizzato
% deriveG(x_) = transpose(jacobian(G(x_(1),x_(2),x_(3),x_(4),x_(5),x_(6)), [x_(1),x_(2),x_(3)]));
% A           = [zeros(3) eye(3);-inv(D(0,0,0,0,0,0))*deriveG(0,0,0,0,0,0) zeros(3)];
% Anumeric    = double(A)
% B           = [zeros(3,1);
%                inv(D(0,0,0,0,0,0))*H]
% Bnumeric    = double(B)
%% matrici sistema non lineare
Anl(x_) = [zeros(3) eye(3);
           zeros(3) -inv(D(x_(1), x_(2), x_(3), x_(4), x_(5), x_(6)))*C(x_(1), x_(2), x_(3), x_(4), x_(5), x_(6))];
Bnl(x_) = [zeros(3,1);
           inv(D(x_(1),x_(2),x_(3),x_(4),x_(5),x_(6)))*H];
Hnl(x_) = [zeros(3,1);
          -inv(D(x_(1), x_(2), x_(3), x_(4), x_(5), x_(6)))*G(x_(1), x_(2), x_(3), x_(4), x_(5), x_(6))];
%% Full nonlinear dynamics for Double Inverted Pendulum system (I guess!)
nonLeqn = Anl(x_(1), x_(2), x_(3), x_(4), x_(5), x_(6))*[x_(1); x_(2); x_(3); x_(4); x_(5); x_(6)] + Hnl(x_(1), x_(2), x_(3), x_(4), x_(5), x_(6)) + Bnl(x_(1), x_(2), x_(3), x_(4), x_(5), x_(6))*sym('u');
%% Linearization of nonlinear system using Jacobian method
Aa(x_)  = jacobian(nonLeqn, [x_(1), x_(2), x_(3), x_(4), x_(5), x_(6)]);
Ass     = double(Aa(0, 0, 0, 0, 0, 0));

Bb(x_)  = jacobian(nonLeqn, sym('u'));
Bss     = double(Bb(0, 0, 0, 0, 0, 0));

%% Progettazione controllore LQR
% Bryson rule: Bryson's method [17] was introduced to adjust Q and R weighting
% matricesas an attempt to overcome the drawbacks of using the trial and error method. 
% According to this rule, Q and R can be calculated by finding the reciprocals
% of squares of the maximum acceptable values of the state and the input control variables. 
Q           = diag([50 50 50 700 700 700]);
R           = 1;
sys         = ss(Ass, Bss, eye(6), 0);
% Co   = ctrb(sys);
% res  = rank(Co); %rank is equal to number of state 6,system is controllable
% Ob   = obsv(sys);
% res2 = rank(Ob); %%rank is equal to number of state 6,system is observable
[K, S, P]   = lqr(sys,Q,R); %p are closed-loop poles system
newsys      = ss((Ass-Bss*K), Bss, eye(6), 0);
%%info = stepinfo(newsys);
%%settling_time = info.SettlingTime; % Tempo di assestamento
% Pole Representation
reale = real(P);
immaginario = imag(P);
% Plot dei poli nel piano complesso (cartesiano)
figure;
plot(reale, immaginario, 'rx', 'MarkerSize', 10); % 'rx' indica crocette rosse
xlabel('Parte Reale');
ylabel('Parte Immaginaria');
title('Diagramma dei Poli (Cartesiano)');
grid on;
axis equal;

%% Simulate the linear system for initial conditions
[y,t,~] = initial(newsys, [0; deg2rad(20); -deg2rad(20); 0; 0; 0], 10);
%[y,t,x]=initial(sysnew,[0;deg2rad(20);deg2rad(20);0;0;0],10);

%%I NEED HERE TO EVALUATE TIME TO SHOW IF SYSTEM IS FASTER OR SLOWER
%%CHANGING Q AND R

% Convert from rad to deg
y(:,2:3)=rad2deg(y(:,2:3));
y(:,5:6)=rad2deg(y(:,5:6));
figure
tl1 = tiledlayout(2,3, 'TileSpacing', 'Compact');
for i=1:1:6
    nexttile
    plot(t,y(:,i)), title('x'+string(i), 'FontSize', 10), 
    grid
end
title(tl1, 'Linear Double-Pendulum System under LQR control')
xlabel(tl1, 'Time')


%% ode45 for nonlinear model
x0     = [0; deg2rad(20); -deg2rad(20); 0; 0; 0]; %initial condition
t_a    = linspace(0, 10,101); %time vector
[t, x] = ode45(@(t, x) nlode(t, x, K, Anl, Bnl, Hnl), t_a, x0);
% Convert from rad to deg
x(:,2:3) = rad2deg(x(:,2:3));
x(:,5:6) = rad2deg(x(:,5:6));
figure
tl2 = tiledlayout(2,3, 'TileSpacing', 'Compact');
for i=1:1:6
    nexttile
    plot(t, x(:,i), 'color', "#D95319"), title('x'+string(i), 'FontSize', 10)
    grid
end
title(tl2, 'Nonlinear Double-Pendulum System under LQR control')
xlabel(tl2,'Time')

and the function is :

%% Nonlinear dynamics of Double Inverted Pendulum on a Cart system
function dxdt = nlode(t, x, K, Anl, Bnl, Hnl)
    % Matrices
    A     = Anl(x(1), x(2), x(3), x(4), x(5), x(6));
    B     = Bnl(x(1), x(2), x(3), x(4), x(5), x(6));
    H     = Hnl(x(1), x(2), x(3), x(4), x(5), x(6));
    
    % LQR controller
    u     = - K*x;  % originally designed based on the Linearized pendulum system
    
    % Nonlinear dynamics
    dxdt_ = A*x + H + B*u;
    dxdt  = double(dxdt_(:));
end

also i want to know why if i change angles and i try higher angles my code became very slow and i cant compute,thanks

[Gianluca Di pietro]

I have to evaluate the error for each state changing Q and R and i have to show that a response is faster or slower changing Q and R

[Sam Chak]

@Gianluca Di pietro,

I'm not entirely clear about the nature of your problem. However, it seems like you're trying to create a loop structure in MATLAB code. Instead of exhaustively analyzing all possible combinations of Q and R, which would be a daunting task due to the infinite number of combinations, why not use a single set of fixed Q and R values that satisfy the desired performance requirements?

If your intention is to demonstrate the impact of adjusting Q and R, I suggest considering the following four possibilities:

1. Significantly increase both Q and R.
2. Significantly increase Q while significantly decreasing R.
3. Significantly decrease Q while significantly increasing R.
4. Significantly decrease both Q and R.

Comparing the results of these four scenarios against the baseline design, where , should provide valuable insights.

[Gianluca Di pietro]

Hello thanks for the answer.I try to be more clear.My aim is to show the effect of changing Q and R on my system.But i have to prove not only with a simple plot..For example i increase Q and decrease R and i see that my sistem became faster..Now i can show that for example considering dominant poles…i saw that increasing Q and decreasing R dominant pole moves to left and system became faster..i also tried to find and compare the settling time from my plots.I found the dominant poles and i found the constant time and then i did 4/5 Times the constant to find the settling time.I wanna know if the method i m following is good to show what we can easily see on plot.So i basically trying to show for example why a system became faster analizying dominant poles and trying to find constant Times from dominant pole…I found that increasing Q an d decreasing R became faster but not for all values for example sometimes looks like increasing Q not always give that faster response.

[Gianluca Di pietro]

I have to create in Matlab a code that could help me to find for each choice of Q and R,the position of poles of my system to show how they change and how they change speed of my system,and I also have to do something to calculate the settling time and constant time system from my states plot to show that High Q makes the system faster and also that stady error should be theorically low comparing to a Low Q.So I need to evaluate stady state error,settling time,constant time and dominant poles to show effect of Q and R.

[Sam Chak]

Hi @Gianluca Di pietro

Could you please provide more details on what you mean by "to find for each choice of Q and R"? What specific criteria or requirements are you looking to satisfy? Additionally, it would be helpful to know how you would like the MATLAB code to approach and handle this task.

[Gianluca Di pietro]

Oh yes i mean that everytime I change value of Q and R in this code I would like to be possible to evaluate settling time and time constant and steady state error of my plots in order to show considering these parameters what Q and R choice changes on my system

[Gianluca Di pietro]

I did some simulations,I started from

0) Q=diag(50) and R=50 then

1) I tried Q=diag(700) and R =100 so I increased both Q and R and I had a system faster and with a higher steady state error than 0)

2) I tried Q=diag(700) and R=1 so i increased Q and decrease R and I had a system faster than 0) and 1) and also with a less steady state error than 0) and 1)

3) Q=diag(1) R =100 I had a system much slower than 0) and 1) and 2) but with smaller steady state error than 0) 1) and 2) …and I expected to be slower but why error was slower with small Q?i expected higher error on final value with low Q

4)Q=diag(10) R=1 so I decrease both Q and R and I had a system faster than 0) and with a higher steady state error than 0)

Answer 1

[Sam Chak]

Hi @Gianluca Di pietro

To observe the state responses, you can manually adjust the design values for the state and input-weighted factors Q and R. Once you have extracted the response characteristics, such as Settling Time and Percent Overshoot, from the response data using the 'stepinfo()' function, you can utilize the 'struct2table()' command to create a table for displaying and comparing these characteristics.

In the example below, I have demonstrated the comparison for and . However, I'm unsure about what you meant by the term "time constant" since it is typically used in reference to first-order systems, which is not the case for your double-inverted pendulum system.

Update: Demonstrated the computation of steady-steady errors and added them into the stepinfo Table.

%% Parameters
m0      = 1.5;
m1      = 0.5;
m2      = 0.75;
L1      = 0.5;
L2      = 0.75;
g       = 9.8;
d1      = m0 + m1 + m2;
d2      = (0.5*m1 + m2)*L1;
d3      = 0.5*m2*L2;
d4      = (1/3*m1 + m2)*L1^2;
d5      = 0.5*m2*L1*L2;
d6      = 1/3*m2*L2^2;
f1      = (0.5*m1 + m2)*L1*g;
f2      = 0.5*m2*L2*g;

%% Symbolic matrices in the nonlinear dynamics
syms Theta0 Theta1 Theta2 DTheta0 DTheta1 DTheta2;
x_      = transpose([Theta0 Theta1 Theta2 DTheta0 DTheta1 DTheta2]);
D(x_)   = [d1,            d2*cos(x_(2)),       d3*cos(x_(3));
           d2*cos(x_(2)), d4,                  d5*cos(x_(2)-x_(3));
           d3*cos(x_(3)), d5*cos(x_(2)-x_(3)), d6];
G(x_)   = [0; 
          -f1*sin(x_(2)); 
          -f2*sin(x_(3))];
C(x_)   = [0, -d2*sin(x_(2))*x_(5),      -d3*sin(x_(3))*x_(6);
           0,  0,                         d5*sin(x_(2)-x_(3))*x_(6);
           0, -d5*sin(x_(2)-x_(3))*x_(5), 0];
H       = [1; 
           0; 
           0];

%% matrici sistema non lineare
Anl(x_) = [zeros(3) eye(3);
           zeros(3) -inv(D(x_(1), x_(2), x_(3), x_(4), x_(5), x_(6)))*C(x_(1), x_(2), x_(3), x_(4), x_(5), x_(6))];
Bnl(x_) = [zeros(3,1);
           inv(D(x_(1),x_(2),x_(3),x_(4),x_(5),x_(6)))*H];
Hnl(x_) = [zeros(3,1);
          -inv(D(x_(1), x_(2), x_(3), x_(4), x_(5), x_(6)))*G(x_(1), x_(2), x_(3), x_(4), x_(5), x_(6))];

%% Full nonlinear dynamics for Double Inverted Pendulum system
nonLeqn = Anl(x_(1), x_(2), x_(3), x_(4), x_(5), x_(6))*[x_(1); x_(2); x_(3); x_(4); x_(5); x_(6)] + Hnl(x_(1), x_(2), x_(3), x_(4), x_(5), x_(6)) + Bnl(x_(1), x_(2), x_(3), x_(4), x_(5), x_(6))*sym('u');

%% Linearization of nonlinear system using Jacobian method
Aa(x_)  = jacobian(nonLeqn, [x_(1), x_(2), x_(3), x_(4), x_(5), x_(6)]);
Ass     = double(Aa(0, 0, 0, 0, 0, 0));
Bb(x_)  = jacobian(nonLeqn, sym('u'));
Bss     = double(Bb(0, 0, 0, 0, 0, 0));

%% Progettazione controllore LQR
Q       = diag([50 50 50 700 700 700]);     % <-- design value input by the User
R       = 1;                                % <-- design value input by the User
sys     = ss(Ass, Bss, eye(6), 0);
[K,S,P] = lqr(sys, Q, R);   % p are closed-loop poles system
newsys  = ss(Ass-Bss*K, Bss, eye(6), 0);

%% Pole Representation
figure;
myLocus = ss(Ass-Bss*K, Bss, [1 0 0 0 0 0], 0);
rlocus(myLocus), grid on, ylim([-15, 15])   % maybe making the graph more meaningful

%% Simulate the linear system for initial conditions
[y, t, ~]   = initial(newsys, [0; deg2rad(10); -deg2rad(10); 0; 0; 0], 5);
y(:,2:3)    = rad2deg(y(:,2:3));    % Convert from rad to deg
y(:,5:6)    = rad2deg(y(:,5:6));

figure
tl1         = tiledlayout(2,3, 'TileSpacing', 'Compact');
for i=1:1:6
    nexttile
    plot(t,y(:,i)), title('x'+string(i), 'FontSize', 10), 
    grid
end
title (tl1, 'Linear Double-Pendulum System under LQR control')
xlabel(tl1, 'Time')

%% Computation of Steady-state Errors (theoretical results)
setpoint= 0;

theta1  = ss(Ass-Bss*K, Bss, [0 1 0 0 0 0], 0); % 2nd state as output
ssOut1  = dcgain(theta1)        % steady-state response of theta 1

ssOut1 = 0

ssErr1  = setpoint - ssOut1     % steady-state error for theta 1

ssErr1 = 0

theta2  = ss(Ass-Bss*K, Bss, [0 0 1 0 0 0], 0); % 3rd state as output
ssOut2  = dcgain(theta2)        % steady-state response of theta 2

ssOut2 = 0

ssErr2  = setpoint - ssOut2     % steady-state error for theta 2

ssErr2 = 0

%% Evaluate Response Characteristics from Response Data
y1                  = y(:,2);   % theta 1
y1final             = y1(end);  % Beware!: not exactly relative to 0°
y1init              = 10;
S1                  = stepinfo(y1, t, y1final, y1init);
S1.SteadyStateError = ssErr1;   % add Steady-state Error into S1 structure

y2                  = y(:,3);   % theta 2
y2final             = y2(end);  % Beware!: not exactly relative to 0°
y2init              = 10;
S2                  = stepinfo(y2, t, y2final, y2init);
S2.SteadyStateError = ssErr2;   % add Steady-state Error into S2 structure

%% Create a stepInfoTable to tabulate the Response Characteristics
stepInfoTable = struct2table([S1, S2]);
stepInfoTable = removevars(stepInfoTable, {'SettlingMin', 'SettlingMax', 'Undershoot', 'Peak', 'PeakTime'});
stepInfoTable.Properties.RowNames = {'theta 1', 'theta 2'};
stepInfoTable

stepInfoTable = 2×5 table
               RiseTime    TransientTime    SettlingTime    Overshoot    SteadyStateError
               ________    _____________    ____________    _________    ________________

    theta 1    0.032158       2.4092           2.5473        180.13             0        
    theta 2           0        2.623           2.6173        97.398             0        

[Gianluca Di pietro]

Yes you are right,thank you so much!I just want to analyze the transient time and the settling time to show how system is becaming faster or slower changing Q and R and prove that also showing how poles moving Also I want if possible to calculate the steady state error respect to my aim of 0º cause usually high Q should provide less deviation from my desired output of 0º..I evaluate my output at settling times and then i calculate the error as the abs(y(settling_time)) since y_des = 0..am I doing right?

[Sam Chak]

Hi @Gianluca Di pietro

The evaluation of steady-state error is typically not performed at the settling time because the term 'steady-state' implies the absence of oscillations. However, the convergence process sometimes takes longer than the simulation time. Therefore, we can use the 'dcgain()' command to theoretically determine the steady-state output. This approach yields results similar to those obtained by applying the Final Value Theorem.

I have also updated the code in my Answer to demonstrate how to compute the steady-state errors for both and . If you find the solution helpful, please consider clicking 'Accept' ✔ on the answer and voting 👍 for it. Your support is greatly appreciated!

[Sam Chak]

Hmm... It's not always zero. The steady-state error of the Cart Position is relatively large, at approximately 0.86. What does this value imply? It means that when a unit step command signal (1 m in the forward direction) is applied to the Cart–Double–Pendulum system, the Cart will eventually stop at a distance of 0.14 m. Therefore, there is a difference of 0.86 m between the desired distance and the actual distance.

%% Parameters
m0      = 1.5;
m1      = 0.5;
m2      = 0.75;
L1      = 0.5;
L2      = 0.75;
g       = 9.8;
d1      = m0 + m1 + m2;
d2      = (0.5*m1 + m2)*L1;
d3      = 0.5*m2*L2;
d4      = (1/3*m1 + m2)*L1^2;
d5      = 0.5*m2*L1*L2;
d6      = 1/3*m2*L2^2;
f1      = (0.5*m1 + m2)*L1*g;
f2      = 0.5*m2*L2*g;

%% Symbolic matrices in the nonlinear dynamics
syms Theta0 Theta1 Theta2 DTheta0 DTheta1 DTheta2;
x_      = transpose([Theta0 Theta1 Theta2 DTheta0 DTheta1 DTheta2]);
D(x_)   = [d1,            d2*cos(x_(2)),       d3*cos(x_(3));
           d2*cos(x_(2)), d4,                  d5*cos(x_(2)-x_(3));
           d3*cos(x_(3)), d5*cos(x_(2)-x_(3)), d6];
G(x_)   = [0; 
          -f1*sin(x_(2)); 
          -f2*sin(x_(3))];
C(x_)   = [0, -d2*sin(x_(2))*x_(5),      -d3*sin(x_(3))*x_(6);
           0,  0,                         d5*sin(x_(2)-x_(3))*x_(6);
           0, -d5*sin(x_(2)-x_(3))*x_(5), 0];
H       = [1; 
           0; 
           0];

%% matrici sistema non lineare
Anl(x_) = [zeros(3) eye(3);
           zeros(3) -inv(D(x_(1), x_(2), x_(3), x_(4), x_(5), x_(6)))*C(x_(1), x_(2), x_(3), x_(4), x_(5), x_(6))];
Bnl(x_) = [zeros(3,1);
           inv(D(x_(1),x_(2),x_(3),x_(4),x_(5),x_(6)))*H];
Hnl(x_) = [zeros(3,1);
          -inv(D(x_(1), x_(2), x_(3), x_(4), x_(5), x_(6)))*G(x_(1), x_(2), x_(3), x_(4), x_(5), x_(6))];

%% Full nonlinear dynamics for Double Inverted Pendulum system
nonLeqn = Anl(x_(1), x_(2), x_(3), x_(4), x_(5), x_(6))*[x_(1); x_(2); x_(3); x_(4); x_(5); x_(6)] + Hnl(x_(1), x_(2), x_(3), x_(4), x_(5), x_(6)) + Bnl(x_(1), x_(2), x_(3), x_(4), x_(5), x_(6))*sym('u');

%% Linearization of nonlinear system using Jacobian method
Aa(x_)  = jacobian(nonLeqn, [x_(1), x_(2), x_(3), x_(4), x_(5), x_(6)]);
Ass     = double(Aa(0, 0, 0, 0, 0, 0));
Bb(x_)  = jacobian(nonLeqn, sym('u'));
Bss     = double(Bb(0, 0, 0, 0, 0, 0));

%% Progettazione controllore LQR
Q       = diag([50 50 50 700 700 700]);     % <-- design value input by the User
R       = 1;                                % <-- design value input by the User
sys     = ss(Ass, Bss, eye(6), 0);
[K,S,P] = lqr(sys, Q, R);   % p are closed-loop poles system
newsys  = ss(Ass-Bss*K, Bss, eye(6), 0);

%% Compute steady-state error for state 1 (x1)
CartPos = ss(Ass-Bss*K, Bss, [1 0 0 0 0 0], 0); % Cart Position as output
ssOut   = dcgain(CartPos)       % steady-state response of Cart Position

ssOut = 0.1414

setpoint= 1;                    % to inject a unit step input signal
step(CartPos), grid on

ssErr   = setpoint - ssOut      % steady-state error for Cart Position

ssErr = 0.8586

[Gianluca Di pietro]

So every steady state error is 0 using the code and I think is right cause my control should be able to take output to 0..maybe when on book says “Q high reduced deviation from desired value” it means that final value will be 0 but maybe oscillations around the desired values are less than in other cases in the sense that I should follow with more precision the desired final state.So to show this is possible the same to use settling time?I mean I decide a tollerance for example 3% around desired output and I choose that it would be my settling time,so for example when values are in this band I have reached the settling time..But of course not all the states changing R and Q are equals maybe some of them are in the 3% others others are in 2% others are closer.So I decided a tolerance to define settling times for all states and for all Q and R and then I see changing Q and R when values in this band are closer than other to final value to see if during time their distance from settling time to “infinity” is closer

[Sam Chak]

Hi @Gianluca Di pietro

However, the steady-state error of the Cart Position in this case can be easily addressed by incorporating a sort of "Pre-amp" at the Reference Input port. This will amplify the commanded signal to a level sufficient for achieving zero steady-state error.

%% Compute steady-state error for state 1 (x1)
Proxy1  = ss(Ass-Bss*K, Bss, [1 0 0 0 0 0], 0);         % Proxy Cart Position model
Preamp  = 1/dcgain(Proxy1)

Preamp = 7.0711

CartPos = ss(Ass-Bss*K, Bss*Preamp, [1 0 0 0 0 0], 0);  % Cart Position as output
ssOut   = dcgain(CartPos)       % steady-state response of Cart Position

ssOut = 1.0000

setpoint= 1;                    % to inject a unit step input signal
step(CartPos), grid on

ssErr   = setpoint - ssOut      % steady-state error for Cart Position

ssErr = 1.1102e-16