How can I reduce oscilltions in control signal in LQI controller and why response is so similar even using different models?

63 views (last 30 days)
Wiktoria
Wiktoria on 28 Oct 2024 at 3:07
Edited: Pavl M. on 28 Oct 2024 at 17:53
I am trying to control a servomechanism that creates pressure by pulling a brake pedal. I have identified the system using System Identification Toolbox. I got the state space matrices, implemented an LQI controller and although the system works, the control signal oscillates when the impulse comes, even when changing reference signal from for example 5 bars to 0.
Although the oscillations in the beginning are small, I want to reduce them as much as possivle, since if the reference signal gets bigger, so do these oscillations. I have tried changing the model I am using and I tried to identify the system again, but surprisingly, this seems to always happen. I am confused why, when I change a model, the behavior is the same. Of course I tried changing the weight matrices Q and R, but even when the oscillations got lower, time of regulation got significantly longer and I want to avoid that too. I want to implement it "step by step", as I want to rewrite it later on to C.
I am completely new to control theory and I wonder what I am doing wrong, so I will be grateful for any help. My current code:
A = [-0.4628, 0.3217, 0.5789, -0.3823, -0.0507;
0.7593, 1.0396, -0.2008, -0.2438, -0.0678;
0.2106, 0.6490, 0.3973, 0.8735, 0.2070;
0.3128, 0.3380, 0.1524, 0.0914, -0.1088;
0.0353, 0.0419, -0.0355, 0.0682, 0.7576];
B = [0.0026;
0.0088;
-0.0036;
0.0490;
-0.1443];
C = [1.1099, 1.3983, 0.8999, -0.0136, 0.0652];
D = 0;
sys = ss(A, B, C, D, 0.01);
step(sys);
[A_size_row, A_size_column] = size(A);
Q = 4*eye(A_size_row+1);
R = 1;
A_aug = [A, zeros(A_size_column, 1); -C, 0];
B_aug = [B; 0];
[K, ~, ~] = lqr(A_aug, B_aug, Q, R);
K_i = K(6);
dt = 0.02;
T = 100; %in miliseconds
time = 0:dt:T;
x = [0;
0;
0;
0;
0];
x_i = 0;
y = [];
u = [];
ref = [];
states = [];
r=0;
for t = time
if t == 50
r = 5;
end
u_t = -K_x * x - K_i * x_i;
x_dot = A * x + B * u_t;
x = x + x_dot * dt;
y_t = C * x;
x_i_dot = r - y_t;
x_i = x_i + x_i_dot * dt;
y = [y, y_t];
u = [u, u_t];
ref = [ref, r];
states = [states, x];
end
figure;
subplot(3,1,1);
plot(time, ref, 'r--', 'LineWidth', 1.5);
hold on;
plot(time, y, 'b', 'LineWidth', 1.5);
xlabel('Time (miliseconds)');
ylabel('Output y');
title('System Output vs Reference Signal');
legend('Reference', 'Output');
grid on;
subplot(3,1,2);
plot(time, u, 'g', 'LineWidth', 1.5);
xlabel('Time (miliseconds)');
ylabel('Control Effort u');
title('Control Effort');
grid on;
subplot(3,1,3);
plot(time, states, 'LineWidth', 1.5);
xlabel('Time (miliseconds)');
ylabel('States');
title('States over Time');
legend('x1', 'x2', 'x3', 'x4', 'x5');
grid on;
  1 Comment
Aquatris
Aquatris on 28 Oct 2024 at 15:20
Edited: Aquatris on 28 Oct 2024 at 15:27
Ran in:
Sometimes you cannot achieve everything you want via a particualar controller structure. So sometimes it is a good idea to play around with how you set your reference. So oneway to deal with it is, since your system is slow enough, you can apply a ramp input that reaches to the end value within 10 seconds instead of step:
A = [-0.4628, 0.3217, 0.5789, -0.3823, -0.0507;
0.7593, 1.0396, -0.2008, -0.2438, -0.0678;
0.2106, 0.6490, 0.3973, 0.8735, 0.2070;
0.3128, 0.3380, 0.1524, 0.0914, -0.1088;
0.0353, 0.0419, -0.0355, 0.0682, 0.7576];
B = [0.0026;
0.0088;
-0.0036;
0.0490;
-0.1443];
C = [1.1099, 1.3983, 0.8999, -0.0136, 0.0652];
D = 0;
sys = ss(A, B, C, D, 0.01);
% step(sys);
[A_size_row, A_size_column] = size(A);
Q = 4*eye(A_size_row+1);
R = 1;
A_aug = [A, zeros(A_size_column, 1); -C, 0];
B_aug = [B; 0];
[K, ~, ~] = lqr(A_aug, B_aug, Q, R);
K_i = K(6);
K_x = K(1:5);
dt = 0.02;
T = 100; %in miliseconds
time = 0:dt:T;
x = [0;
0;
0;
0;
0];
x_i = 0;
y = [];
u = [];
ref = [];
states = [];
r=0;
for t = time
if t >= 50 & t <= 60
r = (t-50)*.5;
elseif t>= 60
r = 5;
else
r = 0;
end
u_t = -K_x * x - K_i * x_i;
x_dot = A * x + B * u_t;
x = x + x_dot * dt;
y_t = C * x;
x_i_dot = r - y_t;
x_i = x_i + x_i_dot * dt;
y = [y, y_t];
u = [u, u_t];
ref = [ref, r];
states = [states, x];
end
figure;
subplot(3,1,1);
plot(time, ref, 'r--', 'LineWidth', 1.5);
hold on;
plot(time, y, 'b', 'LineWidth', 1.5);
xlabel('Time (miliseconds)');
ylabel('Output y');
title('System Output vs Reference Signal');
legend('Reference', 'Output');
grid on;
subplot(3,1,2);
plot(time, u, 'g', 'LineWidth', 1.5);
xlabel('Time (miliseconds)');
ylabel('Control Effort u');
title('Control Effort');
grid on;
ylim([-10 10])
subplot(3,1,3);
plot(time, states, 'LineWidth', 1.5);
xlabel('Time (miliseconds)');
ylabel('States');
title('States over Time');
legend('x1', 'x2', 'x3', 'x4', 'x5');
grid on;

Sign in to comment.

Sign in to answer this question.

Accepted Answer

Pavl M.
Pavl M. on 28 Oct 2024 at 17:49
Edited: Pavl M. on 28 Oct 2024 at 17:53
Ran in:
Good day.
The topic, theme, the agenda of the special regulator control strategy design, verification and implementation for custom actuators and extended pilot control is no doubt very valuable, interesting and I commenced shaping, sharpening, further enhancing it within the process in progress.
The scope at hand actually opens quite valuable large field control area of betterment workouting.
Let's make optimal controller for it,
What are the constraints on control effort(u) energy and speed?
What are the required settling time, overshoot, rise time?
As soon as I complete, contact me for results.
So far solved:
Found:
Whether the system is stabilizable:
In MPPL NCE TCE:
disp('Whether the system is stabilizable:')
isstabilizable(sys)
isstable(sys)
disp('The pair (A,B) must be stabilizable')
ans = 1
ans = 1
The pair (A,B) must be stabilizable.
The system is stable and stabilizable and so optimal Ricatti control exist.
Just to find it by finding the optimal Ki and Kx
1. In your code you actually activated Linear Quadratic Regulator (LQR) lqr function and not Linear Quadratic Integral (LQI) lqi methode.
2. For what to augment the matrices?
3. Try also good: reg = lqg(sys,QXU,QWV,QI,'1dof')
4. Why not to use just PID or fractional PI controller for pressure control?
5. Let me make also Simulink model for it:
It for sure widens and limits the solution space to find the optimal controller for large familty of plants.
Upon working on your code, the osciallation alleviated, just large signals appeared to fix.
So far I've constructed Gaussian Servo 1DOF controller clsys. Hope this to continue to contribute positive value to solving.
Also in this code, which is just a begining of more valuable commercial project I propose to create from it, there is no oscillations in control signal:
I created Simulink structure for it. I think I can do it less oscillatory, more stable. It must needs to get necessary support investment from industry to me, to make it more. Listen, see, reasonable: For more elegant, precise and professional completed results (I can do if you hire / pay me) more substantial funds investment is necessary
clc
clear all
close all
format long
%pkg load control
A = [-0.4628, 0.3217, 0.5789, -0.3823, -0.0507;
0.7593, 1.0396, -0.2008, -0.2438, -0.0678;
0.2106, 0.6490, 0.3973, 0.8735, 0.2070;
0.3128, 0.3380, 0.1524, 0.0914, -0.1088;
0.0353, 0.0419, -0.0355, 0.0682, 0.7576];
B = [0.0026;
0.0088;
-0.0036;
0.0490;
-0.1443];
C = [1.1099, 1.3983, 0.8999, -0.0136, 0.0652];
D = 0;
dt = 0.02;
aug_level = 0;
gain = 2;
sys = ss(A, B, C, D, dt,'InputName','u','OutputName',{'y'});
figure
step(sys);
%isstabilizable(sys)
%isprop(sys)
isstable(sys)
ans = logical
1
disp('The pair (A,B) must be stabilizable')
The pair (A,B) must be stabilizable
[A_size_row, A_size_column] = size(A);
Q1 = gain*eye(A_size_row+aug_level+1);
Q2 = gain*eye(A_size_row+aug_level);
R = 1; %eye(size(A)) requires square of size 1;
if aug_level == 0
A_aug = A;
else
A_aug = [A, zeros(A_size_column, aug_level); -C, zeros(1,aug_level);zeros(aug_level-1,A_size_column+aug_level)]
end
B_aug = [B; zeros(aug_level,1)]
B_aug = 5×1
0.002600000000000 0.008800000000000 -0.003600000000000 0.049000000000000 -0.144300000000000
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
[K1,XG1,L1] = lqi(sys,Q1,R)
K1 = 1×6
1.705594665070935 4.107133993198723 -0.091632065921366 -2.023673740543694 -1.053262715993665 -1.330175460064640
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
XG1 = 6×6
1.097009009585034 2.542373097035204 -0.059251998389178 -1.145298679922787 -0.347999282268025 -0.813564314191786 2.542373097035204 6.097892816111214 -0.125757189583362 -2.776368877554967 -0.843475004002272 -1.933609923156399 -0.059251998389178 -0.125757189583362 0.049314619962354 0.086754761074918 0.030004947734360 0.020372343906181 -1.145298679922787 -2.776368877554967 0.086754761074918 1.386714427976461 0.422100942957140 0.867631799350868 -0.347999282268025 -0.843475004002272 0.030004947734360 0.422100942957140 0.173836006106934 0.265732457405469 -0.813564314191786 -1.933609923156399 0.020372343906181 0.867631799350868 0.265732457405469 1.468714265499774
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
L1 =
0.983451353581629 + 0.000000000000000i 0.910180000501835 + 0.099143266329425i 0.910180000501835 - 0.099143266329425i 0.679803075624528 + 0.000000000000000i -0.377123713773861 + 0.295778081512648i -0.377123713773861 - 0.295778081512648i
%For Continuous system:
%[X2, L2, G2] = care(A_aug,B_aug,Q2,R)
%for discrete system:
[X3,L3,G3] = dare(A_aug,B_aug,Q2,R)
X3 = 5×5
0.597840613293519 1.353153935314747 -0.041830397919194 -0.603858229598839 -0.180856114140849 1.353153935314747 3.264462434344498 -0.083836662547745 -1.485622614971075 -0.444878198713507 -0.041830397919194 -0.083836662547745 0.047950843260205 0.066365519788052 0.023448492376779 -0.603858229598839 -1.485622614971075 0.066365519788052 0.796443982460092 0.239346417709626 -0.180856114140849 -0.444878198713507 0.023448492376779 0.239346417709626 0.117146539048524
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
L3 =
0.909920602241653 + 0.098788758873526i 0.909920602241653 - 0.098788758873526i 0.679793841417242 + 0.000000000000000i -0.377123714887750 + 0.295778083416622i -0.377123714887750 - 0.295778083416622i
G3 = 1×5
0.868875107432436 2.111484330096811 -0.058292605736457 -1.107623528381077 -0.768786874572413
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
[K,XG,L] = dlqr(A_aug, B_aug, Q2, R)
K = 1×5
0.868875107432408 2.111484330096744 -0.058292605736450 -1.107623528381036 -0.768786874572399
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
XG = 5×5
0.597840613293508 1.353153935314717 -0.041830397919192 -0.603858229598824 -0.180856114140844 1.353153935314717 3.264462434344424 -0.083836662547739 -1.485622614971036 -0.444878198713493 -0.041830397919192 -0.083836662547739 0.047950843260205 0.066365519788049 0.023448492376778 -0.603858229598824 -1.485622614971036 0.066365519788049 0.796443982460071 0.239346417709619 -0.180856114140844 -0.444878198713493 0.023448492376778 0.239346417709619 0.117146539048522
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
L =
0.909920602241651 + 0.098788758873526i 0.909920602241651 - 0.098788758873526i 0.679793841417243 + 0.000000000000000i -0.377123714887750 + 0.295778083416623i -0.377123714887750 - 0.295778083416623i
%Optional!, not final, ToDo: in progress
K_i = -min(X3)
K_i = 1×5
0.603858229598839 1.485622614971075 0.083836662547745 1.485622614971075 0.444878198713507
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
K_x = -X3
K_x = 5×5
-0.597840613293519 -1.353153935314747 0.041830397919194 0.603858229598839 0.180856114140849 -1.353153935314747 -3.264462434344498 0.083836662547745 1.485622614971075 0.444878198713507 0.041830397919194 0.083836662547745 -0.047950843260205 -0.066365519788052 -0.023448492376779 0.603858229598839 1.485622614971075 -0.066365519788052 -0.796443982460092 -0.239346417709626 0.180856114140849 0.444878198713507 -0.023448492376779 -0.239346417709626 -0.117146539048524
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
T = 100; %in miliseconds
time = 0:dt:T;
x = [0;0;0;0;0];
x_i = x;
y = [];
u = [];
ref = [];
states = [];
r=0;
for t = time
if t == 50
r = 5;
end
u_t = -K_x*x - K_i*x_i;
u_t(abs(u_t)>abs(max(C))) = max(C);
x(abs(x)>abs(max(r))) = r;
x_dot = A * x + B' * u_t;
x = x + x_dot * dt;
y_t = C * x;
x_i_dot = r - y_t;
x_i = x_i + x_i_dot * dt;
y = [y, y_t];
u = [u, u_t];
ref = [ref, r];
states = [states, x];
end
figure;
subplot(3,1,1);
plot(time, ref, 'r--', 'LineWidth', 1.5);
hold on;
plot(time, y, 'b', 'LineWidth', 1.5);
xlabel('Time (miliseconds)');
ylabel('Output y');
title('System Output vs Reference Signal');
legend('Reference', 'Output');
grid on;
subplot(3,1,2);
plot(time, u, 'g', 'LineWidth', 1.5);
xlabel('Time (miliseconds)');
ylabel('Control Effort u');
title('Control Effort');
grid on;
subplot(3,1,3);
plot(time, states, 'LineWidth', 1.5);
xlabel('Time (miliseconds)');
ylabel('States');
title('States over Time');
legend('x1', 'x2', 'x3', 'x4', 'x5');
grid on;
%N = eye(size(A))
%[K,S,e] = lqry(sys,Q,R,N)
%% More robust gaussian servo controller:
nx = 5; %Number of states
ny = 5; %Number of outputs
%Disturbance parameters, covariance matrix, quite small as in real world:
Qn = [0.04 0.02 0 0 0; 0.02 0.01 0 0 0; 0 0 0.01 0 0; 0 0 0 0 0; 0 0 0 0 0]
Qn = 5×5
0.040000000000000 0.020000000000000 0 0 0 0.020000000000000 0.010000000000000 0 0 0 0 0 0.010000000000000 0 0 0 0 0 0 0 0 0 0 0 0
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
Rn = 0.01;
R = sqrt(Qn); %[0.01 0 0 0 0; 0 0.02 0 0; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0]
QXU = blkdiag(0.1*eye(nx),R)
QXU = 10×10
0.100000000000000 0 0 0 0 0 0 0 0 0 0 0.100000000000000 0 0 0 0 0 0 0 0 0 0 0.100000000000000 0 0 0 0 0 0 0 0 0 0 0.100000000000000 0 0 0 0 0 0 0 0 0 0 0.100000000000000 0 0 0 0 0 0 0 0 0 0 0.200000000000000 0.141421356237310 0 0 0 0 0 0 0 0 0.141421356237310 0.100000000000000 0 0 0 0 0 0 0 0 0 0 0.100000000000000 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
QWV = blkdiag(Qn,Rn)
QWV = 6×6
0.040000000000000 0.020000000000000 0 0 0 0 0.020000000000000 0.010000000000000 0 0 0 0 0 0 0.010000000000000 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.010000000000000
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
QI = 1 ; %eye(ny)
%Servo-controller design:
KLQG1 = lqg(sys,QXU(1:A_size_row+1,1:A_size_column+1),QWV,QI,'1dof')
KLQG1 = A = x1_e x2_e x3_e x4_e x5_e xi1 x1_e -0.3974 0.399 0.6359 -0.378 -0.04456 0.00563 x2_e 0.1016 0.1935 -0.7207 -0.2184 -0.09865 0.01906 x3_e -0.1348 0.221 0.1118 0.8706 0.1835 -0.007796 x4_e -0.04068 -0.2048 -0.06024 0.1925 -0.08624 0.1061 x5_e 0.2555 0.6064 -0.07474 -0.2194 0.643 -0.3125 xi1 0 0 0 0 0 1 B = e1 x1_e 0.06304 x2_e -0.5787 x3_e -0.3169 x4_e -0.2414 x5_e -0.02868 xi1 0.02 C = x1_e x2_e x3_e x4_e x5_e xi1 u -1.747 -4.19 0.09304 1.996 0.7815 2.166 D = e1 u 0 Input groups: Name Channels Error 1 Output groups: Name Channels Controls 1 Sample time: 0.02 seconds Discrete-time state-space model.
KLQG2 = lqg(sys,QXU(1:A_size_row+1,1:A_size_column+1),QWV,QI,'2dof')
KLQG2 = A = x1_e x2_e x3_e x4_e x5_e xi1 x1_e -0.3974 0.399 0.6359 -0.378 -0.04456 0.00563 x2_e 0.1016 0.1935 -0.7207 -0.2184 -0.09865 0.01906 x3_e -0.1348 0.221 0.1118 0.8706 0.1835 -0.007796 x4_e -0.04068 -0.2048 -0.06024 0.1925 -0.08624 0.1061 x5_e 0.2555 0.6064 -0.07474 -0.2194 0.643 -0.3125 xi1 0 0 0 0 0 1 B = r1 y x1_e 0 -0.06304 x2_e 0 0.5787 x3_e 0 0.3169 x4_e 0 0.2414 x5_e 0 0.02868 xi1 0.02 -0.02 C = x1_e x2_e x3_e x4_e x5_e xi1 u -1.747 -4.19 0.09304 1.996 0.7815 2.166 D = r1 y u 0 0 Input groups: Name Channels Setpoint 1 Measurement 2 Output groups: Name Channels Controls 1 Sample time: 0.02 seconds Discrete-time state-space model.
%Plot outcome - control product response:
figure
step(KLQG1)
title('1DOF Servo Controller designed')
%Constructing closed-loop system after designing:
%Ther are many similar questions this answer will fit(answer) to:
%https://www.mathworks.com/matlabcentral/answers/484846-how-to-use-the-generated-lqg-controller
C=KLQG1;
P = sys;
set(C,'InputName','e')
set(C,'OutputName','u')
set(P,'OutputName','y')
set(P,'InputName','u')
Sum_Error = sumblk('e = r - y');
clsys = tf(connect(C,P,Sum_Error,'r','y'))
clsys = From input "r" to output "y": 0.003423 z^9 + 0.00423 z^8 - 0.01212 z^7 - 0.03851 z^6 + 0.1017 z^5 - 0.06887 z^4 - 0.002825 z^3 + 0.01118 z^2 + 0.006176 z - 0.004271 ------------------------------------------------------------------------------------------------------------------------------------------------ z^11 - 3.567 z^10 + 4.632 z^9 - 2.505 z^8 + 0.5232 z^7 - 0.2944 z^6 + 0.2939 z^5 - 0.1068 z^4 + 0.04219 z^3 - 0.02263 z^2 + 0.005732 z - 0.00111 Sample time: 0.02 seconds Discrete-time transfer function.
clsys_m = minreal(clsys)
clsys_m = From input "r" to output "y": 0.003423 z^9 + 0.00423 z^8 - 0.01212 z^7 - 0.03851 z^6 + 0.1017 z^5 - 0.06887 z^4 - 0.002825 z^3 + 0.01118 z^2 + 0.006176 z - 0.004271 ------------------------------------------------------------------------------------------------------------------------------------------------ z^11 - 3.567 z^10 + 4.632 z^9 - 2.505 z^8 + 0.5232 z^7 - 0.2944 z^6 + 0.2939 z^5 - 0.1068 z^4 + 0.04219 z^3 - 0.02263 z^2 + 0.005732 z - 0.00111 Sample time: 0.02 seconds Discrete-time transfer function.
figure
step(clsys_m)
title('Constructed closed loop system response')
x = [0;0;0;0;0];
x_i = x;
y = [];
u = [];
ref = [];
states = [];
r=0;
%Todo: simulate it in loop:
%noise = 0.001;
%for t = time
% if t == 50
% r = 5;
% end
% u_t = ? + noise;
% x_dot = clsys.A * x + clsys.B * u_t + noise;
% x = x + x_dot * dt;
% y_t = KLQG1.C * x + KLQG1.D*u_;
% x_i_dot = r - y_t;
% x_i = x_i + x_i_dot * dt;
% y = [y, y_t];
% u = [u, u_t];
% ref = [ref, r];
% states = [states, x];
%end
%figure;
%subplot(3,1,1);
%plot(time, ref, 'r--', 'LineWidth', 1.5);
%hold on;
%plot(time, y, 'b', 'LineWidth', 1.5);
%xlabel('Time (miliseconds)');
%ylabel('Output y');
%title('System Output vs Reference Signal');
%legend('Reference', 'Output');
%grid on;
%subplot(3,1,2);
%plot(time, u, 'g', 'LineWidth', 1.5);
%xlabel('Time (miliseconds)');
%ylabel('Control Effort u');
%title('Control Effort');
%grid on;
%subplot(3,1,3);
%plot(time, states, 'LineWidth', 1.5);
%xlabel('Time (miliseconds)');
%ylabel('States');
%title('States over Time');
%legend('x1', 'x2', 'x3', 'x4', 'x5');
%grid on;
%Option 2, construct positive semidefinite matrices:
%Nn = 1;
%[kest,L,P2] = kalman(sys,Qn,Rn,Nn);
%Connect the and sys to form the working loop, it should bring more performant control:
%K4 = lqi(sys,Q,R)
%servocontroller = lqgtrack(kest, K4);
%C=servocontroller;
%set(C,'InputName','e')
%set(C,'OutputName','u')
%set(P,'OutputName','y')
%set(P,'InputName','u')
%Sum_Error = sumblk('e = r - y');
%clsys2 = minreal(tf(connect(C,P,Sum_Error,'r','y')))
%figure
%step(clsys2)
%title('Constructed closed loop system response 2')

More Answers (0)

Categories

Find more on General Applications in Help Center and File Exchange

Tags

Products


Release

R2022b

Community Treasure Hunt

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

Start Hunting!