Portfolio
  • Projects
    • Machine Learning
    • Control Systems
    • Robotics
    • Machine Vision
  • About
  • Work Experience
  • Resources
    • My Recommendations
    • Modern Robotics
    • Robot Modelling and Control
    • Python for Data Science
  • Contact
    • LinkedIn
    • GitHub
    • Stack Overflow
    • Gmail
    • Curriculum Vitae
  1. Projects
  2. Control System
  3. Control Design
  • Home
  • Projects
    • Robotics
      • Inverse Kinematics
      • Trajectory Optimization
      • Dynamics
      • Control
      • RRT
    • Control System
      • Control Design
      • LTI Modeling
      • Bode Estimation
    • Machine Learning
      • Probability
      • Regression
      • Clustering
      • Classification
      • Anamoly
    • Machine Vision
      • Fire Detection
  • About Me
  • Tools
    • Git
    • Quarto
  • Resources
    • My Recommendations
    • Modern Robotics
    • Robot Modelling and Control
    • Python for Data Science

Full State Feedback Control Design

Control
Optimal LQR control for drone landing using MATLAB, emphasizing strategic pole placement for precise descent

ME5554/AOE5754/ECE5754 Applied Linear Systems - Project II

clear;
set(0,'defaultlinelinewidth',2);
format short;
sympref('FloatingPointOutput',true);

CASE 1. Prepare the Open-Loop Plant

The state-space LTI object we developed in the Project I for the open-loop plant should have been saved in a .mat file.


Note: Before we proceed with the implementation, one should be responsible for correcting any errors they have in their open-loop plant! In this case, the plant is provided in my GitHub. One should NOT regenerate the open-loop LTI state-space plant in this LiveScript!!! The open-loop plant MUST be loaded from a .mat file (created in ME5554/AOE5754/ECE5754 Applied Linear Systems - Project I - Linearized LTI Modeling Feedback Control) to avoid discrepancies.

load("MP_nishant_ss_model.mat")

After loading the .mat file, display the open-loop numerical state-space model. This model MUST display appropriately chosen statenames, inputnames (ok to use defaults here), and outputnames.

OL_sys
OL_sys =
 
  A = 
                   X Position    X Velocity    Y Position    Y Velocity  Angular Posi  Angular Velo
   X Position               0             1             0             0             0             0
   X Velocity               0            -1             0             0         -9.81             0
   Y Position               0             0             0             1             0             0
   Y Velocity               0             0             0            -3             0             0
   Angular Posi             0             0             0             0             0             1
   Angular Velo             0             0             0             0             0          -1.6
 
  B = 
                 Varitional F  Varitional F
   X Position               0             0
   X Velocity               0             0
   Y Position               0             0
   Y Velocity             0.1           0.1
   Angular Posi             0             0
   Angular Velo        0.1302       -0.0898
 
  C = 
                   X Position    X Velocity    Y Position    Y Velocity  Angular Posi  Angular Velo
   X Position               1             0             0             0             0             0
   Y Position               0             0             1             0             0             0
   Angular Posi             0             0             0             0             1             0
 
  D = 
                 Varitional F  Varitional F
   X Position               0             0
   Y Position               0             0
   Angular Posi             0             0
 
Continuous-time state-space model.
Model Properties

Next, extract the Aol, Bol, Col, and Dol matrices, but do not display these results (because we already displayed them in the previous step).

Aol = OL_sys.A;
Bol = OL_sys.B;
Col = OL_sys.C;
Dol = OL_sys.D;

Compute the open-loop eigenvalues (ol_poles), and the system dimensions (N, M, P). Assign these quantities to variables for use later in the LiveScript, but do not display them here.

ol_poles = eig(Aol)
ol_poles = 6x1    
         0
   -1.0000
         0
   -3.0000
         0
   -1.6000
N = size(Aol,2)
N = 6
M = size(Bol,2)
M = 2
P = size(Col,1)
P = 3

CASE 2. Simulate Closed-Loop Control

We are given the following set of desired closed-loop poles which are known to result in 5% settling times (transient times) that are less than 3 seconds without saturating any of the control signals.

cl_poles = [
    -29 + 15i
    -29 - 15i
    -32 + 0i
    -1.6 + 0.76i
    -1.6 - 0.76i
    -1.2 + 0i];

Use the PLACE function to compute the full state feedback control gain matrix, and display our result.

G = place(Aol,Bol,cl_poles)
G = 2x6    
1.0e+03 *

   -1.1332   -0.8551    0.0362    0.1151    4.2932    0.2707
    2.4067    1.7537    0.4483    0.1343   -8.6992   -0.3186

Assume that the actual plant starts with the following initial conditions (NOTE: These are NEW IC’s!!!):

x0 = -2;            % meters
x0dot = 0.5;        % meters/sec
y0 = 1;             % meters
y0dot = 0;          % meters/sec
theta0 = -0.5;      % radians
theta0dot = 0.1;    % radians/sec

Assume that we have access to the full state vector, then construct a closed-loop state-space LTI object using the feedback gains from above.

Acl = Aol - Bol*G
Acl = 6x6    
1.0e+03 *

         0    0.0010         0         0         0         0
         0   -0.0010         0         0   -0.0098         0
         0         0         0    0.0010         0         0
   -0.1273   -0.0899   -0.0484   -0.0279    0.4406    0.0048
         0         0         0         0         0    0.0010
    0.3637    0.2688    0.0355   -0.0029   -1.3402   -0.0655
Ccl = Col - Dol*G
Ccl = 3x6    
     1     0     0     0     0     0
     0     0     1     0     0     0
     0     0     0     0     1     0

We must augment the output equation to output the original outputs as well as the control signals.

Ccl = [Ccl; -G]
Ccl = 5x6    
1.0e+03 *

    0.0010         0         0         0         0         0
         0         0    0.0010         0         0         0
         0         0         0         0    0.0010         0
    1.1332    0.8551   -0.0362   -0.1151   -4.2932   -0.2707
   -2.4067   -1.7537   -0.4483   -0.1343    8.6992    0.3186
CL_sys = ss(Acl,[],Ccl,[]);

Simulate the initial value response using INITIAL. Choose a final simulation time of 4 seconds for this simulation.

t_final = 4; % seconds
xinit = [x0;x0dot;y0;y0dot;theta0;theta0dot];
[y_ivr,t_ivr] = initial(CL_sys,xinit,t_final);

Estimate the 5% settling times for each output.

S = lsiminfo(y_ivr,t_ivr,"SettlingTimeThreshold",0.05);
S(1:3).TransientTime % Settling times in seconds
ans = 2.4134
ans = 2.6797
ans = 0.9707

Generate a plot like what was required for the Midterm Project with the three outputs on the left-hand side (separate axes for each output), and the control signals on the right-hand side of the figure (separate axes for each control signal). Annotate each plot with some kind of marker for the 5% settling times, the min and max allowable control signals, grid lines, axis labels, and units.

f = figure();

u1_min = -2039;
u1_max = 1959;
u2_min = -2057.1;
u2_max = 1940.9;
outputNames = {'x', 'y', '\theta', 'u1','u2'};
outputUnits = {'m', 'm', 'rad','N', 'N'};
inputNames = {'u1', 'u2'};
inputUnits = {'N', 'N'};

for i = 1:3
    subplot(3, 2, -1+2*i); % odd numbers for left side
  
    plot(t_ivr, y_ivr(:, i));
    
    index = find((t_ivr-S(i).TransientTime)>0,1);
    if index<=size(t_ivr,1)
        xline(t_ivr(index),'--r')
    end
    clear index

    title(sprintf('%s', outputNames{i}));
    xlabel('Time (s)');
    ylabel(sprintf('%s', outputUnits{i}));
    grid on
end

for i = 1:2
    subplot(2, 2, i*2); % even numbers for right side
    
    plot(t_ivr, y_ivr(:, i+3)); % 4 and 5 for u1 and u2
    ylim([-2500,2500])
    if i == 1
        yline(double(u1_max),'--b')
        yline(double(u1_min),'--b')
    end
    if i == 2
        yline(double(u2_max),'--b')
        yline(double(u2_min),'--b')
    end

    title(sprintf('%s', inputNames{(i)}));
    xlabel('Time (s)');
    ylabel(sprintf('%s', inputUnits{(i)}));
    grid on
    
end
% Adjust the overall plot title
sgtitle('Initial Response for Outputs (left) & Inputs (right)');
set(f,"Position", [0 0 800 800])

In addition to the normal five subplot figure, we must also include the animation plot in the LiveScript as documentation of the performance.

als_f23_animate( t_ivr, y_ivr(:,1), y_ivr(:,2), y_ivr(:,3));

CASE 3. Evaluate Observability

First show that this system is completely observable using all three outputs.

Ob = obsv(Aol,Col);
rank(Ob)
ans = 6
N
N = 6

Rank of observability matrix is equal to number of states.

Then we determine whether this system is completely observable using any subset of outputs and justify our answers.

For individual outputs

rank(obsv(Aol,Col(1,:)))
ans = 4
rank(obsv(Aol,Col(2,:)))
ans = 2
rank(obsv(Aol,Col(3,:)))
ans = 2

None of the individual outputs obs matrices are full rank

rank(obsv(Aol,Col([1 2],:)))
ans = 6
rank(obsv(Aol,Col([1 3],:)))
ans = 4
rank(obsv(Aol,Col([2 3],:)))
ans = 4

We can see that the combination of ouput 1 & 2 (x & y), has a obs matrix which is full rank meaning the system is obervable using just there 2 outputs.

CASE 4. Design a Full-Order Luenberger Observer

In this CASE, we will compute the gains for a state estimator using all three outputs. We must first choose appropriate pole locations for the state estimator (observer), and then compute and display the associated feedback gain matrix using the PLACE function.

alpha = 1.5;
obs_poles = alpha*cl_poles;% alpha > 1 is faster observer than controller
K = place(Aol.',Col.',obs_poles).'
K = 6x3    
   47.0621   11.3921  -12.6548
   59.1867  -53.1150   13.2504
   -6.5959   78.7139  -41.1033
   63.5968    2.2363   -0.7294
   -5.7511   40.1340   10.2240
    1.9765   -0.8166    6.4001

To document our design, we must generate a plot that includes the open-loop poles from CASE 1, the closed-loop poles from CASE 2, and observer poles from this CASE. Use common markers for each set of poles (i.e. all open-loop poles must use the same marker, etc.) so that it is easy to distinguish the poles on the figure. Include a legend that only has three entries in it (open-loop, closed-loop, observer), corresponding to the three sets of poles. We must use the SGRID function (read the Matlab help!!) to generate appropriate grid lines. We will also annotate our plot with axis labels. Making sure our legend does not cover up any of the poles!

f = figure;

wmax = round(max(abs([cl_poles;ol_poles;obs_poles])));
sgrid(0:0.1:0.9,0:wmax/5:wmax*1.5,"new")
h1 = plot(real(ol_poles),imag(ol_poles),'bx','markersize',10,'linewidth',2);
hold on
h2 = plot(real(cl_poles),imag(cl_poles),'gs','markersize',10,'linewidth',2);
h3 = plot(real(obs_poles),imag(obs_poles),'rd','markersize',10,'linewidth',2);
xlabel('Real (s)')
ylabel('Imag (s)')
xline(0,'m','linewidth',2)
yline(0,'m','linewidth',2)
legend([h1,h2,h3],'open loop','closed loop','estimator');
set(f,"Position", [0 0 600 600])
set(gca,'yaxislocation','right')
axis equal

CASE 5. Simulate Output Feedback Control

For this CASE, we must simulate the complete output feedback control system, which will include the state estimator and full state feedback using the estimated states.

Using the equations for obersever followed by controller we can find, (L14/S29)

\[ \left\lbrack \begin{array}{c} \dot{x} \\ \dot{\hat{x} } \end{array}\right\rbrack =\left\lbrack \begin{array}{cc} A & -B\;G\\ K\;C & \hat{A} -K\hat{C} -\hat{B\;} G-K\left(D-\hat{D} \right)G \end{array}\right\rbrack \left\lbrack \begin{array}{c} x\\ \hat{x} \end{array}\right\rbrack \]

\[ \left\lbrack \begin{array}{c} y\\ \hat{y} \end{array}\right\rbrack =\left\lbrack \begin{array}{cc} C & -D\;G\\ 0 & \hat{C} -\hat{D} \;G \end{array}\right\rbrack \left\lbrack \begin{array}{c} x\\ \hat{x} \end{array}\right\rbrack \]

Assumuing,

\[ \hat{A} =A\;\;\;\;\;\;\;\;\;\;\;\;\hat{B} =B\;\;\;\;\;\;\;\;\;\;\;\;\hat{C} =C\;\;\;\;\;\;\;\;\;\;\;\;\hat{D} =D \]

We can simplify as,

\[ \left\lbrack \begin{array}{c} \dot{x} \\ \dot{\hat{x} } \end{array}\right\rbrack =\left\lbrack \begin{array}{cc} A & -B\;G\\ K\;C & A-K\;C-B\;G \end{array}\right\rbrack \left\lbrack \begin{array}{c} x\\ \hat{x} \end{array}\right\rbrack \]

\[ \left\lbrack \begin{array}{c} y\\ \hat{y} \end{array}\right\rbrack =\left\lbrack \begin{array}{cc} C & -D\;G\\ 0 & C-D\;G \end{array}\right\rbrack \left\lbrack \begin{array}{c} x\\ \hat{x} \end{array}\right\rbrack \]

Augmenting to get back u to plot,

\[ \begin{bmatrix} y \\ \hat{y} \\ u \end{bmatrix} = \begin{bmatrix} C & -DG \\ 0 & C-DG \\ 0 & -G \end{bmatrix} \begin{bmatrix} x \\ \hat{x} \end{bmatrix} \]

In order to insure that the estimated states \(\hat{\mathbf{x}}\) (and outputs \(\hat{\mathbf{y}}\)) have converged before control is enabled, we will first need to generate an augmented LTI closed-loop output-feedback control model with the state feedback gain \(\mathbf{G}\) set to zero (it is easier to just use \(\left(0*\mathbf{G}\right)\) in the matrices). This first LTI object will be used for the first 200 milliseconds of the response. We must use the same initial states on the plant as defined in CASE 2 above, but the initial observer states must be zero. Our closed-loop state-space output equation must also be augmented to output the following three vectors: \(\mathbf{y}\), \(\hat{\mathbf{y}}\), and \(\mathbf{u}\).

xinit_OFC = [xinit;zeros(size(xinit))];

Aofc_G0 = [Aol, -Bol*G*0;
    K*Col, Aol-K*Col-Bol*G*0];
Cofc_G0 = [Col, -Dol*G*0;
    zeros(size(Col)), Col-Dol*G*0;
    zeros(size(G)), -G*0;
    eye(6,6),zeros(6,6);
    zeros(6,6),eye(6,6)];
OFC_sys_G0 = ss(Aofc_G0,[],Cofc_G0,[])
OFC_sys_G0 =
 
  A = 
             x1       x2       x3       x4       x5       x6       x7       x8       x9      x10      x11      x12
   x1         0        1        0        0        0        0        0        0        0        0        0        0
   x2         0       -1        0        0    -9.81        0        0        0        0        0        0        0
   x3         0        0        0        1        0        0        0        0        0        0        0        0
   x4         0        0        0       -3        0        0        0        0        0        0        0        0
   x5         0        0        0        0        0        1        0        0        0        0        0        0
   x6         0        0        0        0        0     -1.6        0        0        0        0        0        0
   x7     47.06        0    11.39        0   -12.65        0   -47.06        1   -11.39        0    12.65        0
   x8     59.19        0   -53.12        0    13.25        0   -59.19       -1    53.12        0   -23.06        0
   x9    -6.596        0    78.71        0    -41.1        0    6.596        0   -78.71        1     41.1        0
   x10     63.6        0    2.236        0  -0.7294        0    -63.6        0   -2.236       -3   0.7294        0
   x11   -5.751        0    40.13        0    10.22        0    5.751        0   -40.13        0   -10.22        1
   x12    1.976        0  -0.8166        0      6.4        0   -1.976        0   0.8166        0     -6.4     -1.6
 
  B = 
     Empty matrix: 12-by-0
 
  C = 
         x1   x2   x3   x4   x5   x6   x7   x8   x9  x10  x11  x12
   y1     1    0    0    0    0    0    0    0    0    0    0    0
   y2     0    0    1    0    0    0    0    0    0    0    0    0
   y3     0    0    0    0    1    0    0    0    0    0    0    0
   y4     0    0    0    0    0    0    1    0    0    0    0    0
   y5     0    0    0    0    0    0    0    0    1    0    0    0
   y6     0    0    0    0    0    0    0    0    0    0    1    0
   y7     0    0    0    0    0    0    0    0    0    0    0    0
   y8     0    0    0    0    0    0    0    0    0    0    0    0
   y9     1    0    0    0    0    0    0    0    0    0    0    0
   y10    0    1    0    0    0    0    0    0    0    0    0    0
   y11    0    0    1    0    0    0    0    0    0    0    0    0
   y12    0    0    0    1    0    0    0    0    0    0    0    0
   y13    0    0    0    0    1    0    0    0    0    0    0    0
   y14    0    0    0    0    0    1    0    0    0    0    0    0
   y15    0    0    0    0    0    0    1    0    0    0    0    0
   y16    0    0    0    0    0    0    0    1    0    0    0    0
   y17    0    0    0    0    0    0    0    0    1    0    0    0
   y18    0    0    0    0    0    0    0    0    0    1    0    0
   y19    0    0    0    0    0    0    0    0    0    0    1    0
   y20    0    0    0    0    0    0    0    0    0    0    0    1
 
  D = 
     Empty matrix: 20-by-0
 
Continuous-time state-space model.
Model Properties

Next, we must generate a second augmented LTI closed-loop output-feedback control model that is identical to the first model except that the state-feedback control gains are not zero, i.e. the controller is enabled. This second LTI object will be used for the remainder of the simulation.

Aofc = [Aol, -Bol*G;
    K*Col, Aol-K*Col-Bol*G];
Cofc = [Col, -Dol*G;
    zeros(size(Col)), Col-Dol*G;
    zeros(size(G)), -G];
OFC_sys = ss(Aofc,[],Cofc,[])
OFC_sys =
 
  A = 
             x1       x2       x3       x4       x5       x6       x7       x8       x9      x10      x11      x12
   x1         0        1        0        0        0        0        0        0        0        0        0        0
   x2         0       -1        0        0    -9.81        0        0        0        0        0        0        0
   x3         0        0        0        1        0        0        0        0        0        0        0        0
   x4         0        0        0       -3        0        0   -127.3   -89.85   -48.45   -24.94    440.6    4.792
   x5         0        0        0        0        0        1        0        0        0        0        0        0
   x6         0        0        0        0        0     -1.6    363.7    268.8    35.54   -2.918    -1340   -63.86
   x7     47.06        0    11.39        0   -12.65        0   -47.06        1   -11.39        0    12.65        0
   x8     59.19        0   -53.12        0    13.25        0   -59.19       -1    53.12        0   -23.06        0
   x9    -6.596        0    78.71        0    -41.1        0    6.596        0   -78.71        1     41.1        0
   x10     63.6        0    2.236        0  -0.7294        0   -190.9   -89.85   -50.68   -27.94    441.3    4.792
   x11   -5.751        0    40.13        0    10.22        0    5.751        0   -40.13        0   -10.22        1
   x12    1.976        0  -0.8166        0      6.4        0    361.7    268.8    36.36   -2.918    -1347   -65.46
 
  B = 
     Empty matrix: 12-by-0
 
  C = 
           x1      x2      x3      x4      x5      x6      x7      x8      x9     x10     x11     x12
   y1       1       0       0       0       0       0       0       0       0       0       0       0
   y2       0       0       1       0       0       0       0       0       0       0       0       0
   y3       0       0       0       0       1       0       0       0       0       0       0       0
   y4       0       0       0       0       0       0       1       0       0       0       0       0
   y5       0       0       0       0       0       0       0       0       1       0       0       0
   y6       0       0       0       0       0       0       0       0       0       0       1       0
   y7       0       0       0       0       0       0    1133   855.1   -36.2  -115.1   -4293  -270.7
   y8       0       0       0       0       0       0   -2407   -1754  -448.3  -134.3    8699   318.6
 
  D = 
     Empty matrix: 8-by-0
 
Continuous-time state-space model.
Model Properties

We are now ready to simulate the closed-loop response, but it must happen in two separate simulation stages. In the first stage, we will use the INITIAL function to simulate for 200 milliseconds. At the end of this first simulation stage, we must capture the final state vector and use it as the initial state vector for the second simulation stage. The total simulation duration should be 4 seconds like above.

t_break = 200e-03; % seconds
t_final = 4; % seconds
dt = 0.01;
t_1 = (0:dt:t_break);
t_2 = (t_break+dt):dt:t_final;
[y_1,t_1] = initial(OFC_sys_G0,xinit_OFC,t_1);
xinit_OFC_break = y_1(size(y_1,1),9:20).';
[y_2,t_2] = initial(OFC_sys, xinit_OFC_break,t_2);
y_ivr = [y_1(:,1:8);y_2(:,1:8)];
t_ivr = [t_1;t_2];
S = lsiminfo(y_ivr,t_ivr,"SettlingTimeThreshold",0.05);
S(1:3).TransientTime
ans = 2.5721
ans = 1.6099
ans = 2.1822

Generate the standard annotated plot (three output subplots on the left and two control signal subplots on the right), with appropriate annotations; but we must also plot the estimated output signals on top of the associated actual output signals. It must be clear from our plots that the estimated outputs all start from zero and converge to the actual outputs within 200 milliseconds. One must NOT try to get the observer states to converge too fast because it will be difficult to see the convergence on the plots. If we have done this step correctly, we should see 5% settling times that are comparable to the original design from CASE 2. The 5% settling times on the outputs must all be below 3 seconds, and we must not saturate either of the control signals.

f = figure();
subplot(3,2,1)
plot(t_ivr,y_ivr(:,1),'r')
hold on
plot(t_ivr,y_ivr(:,4),'--b')
xline(S(1).TransientTime,'--r')
xlim([0,t_final])
legend({'${x}$','$\hat{x}$'},'Interpreter','latex','Location','best')
ylabel('m')
grid on

subplot(3,2,3)
plot(t_ivr,y_ivr(:,2),'g')
hold on
plot(t_ivr,y_ivr(:,5),'--black')
xlim([0,t_final])
xline(S(2).TransientTime,'--r')
legend({'${y}$','$\hat{y}$'},'Interpreter','latex','Location','best')
ylabel('m')
grid on

subplot(3,2,5)
plot(t_ivr,y_ivr(:,3),'m')
hold on
plot(t_ivr,y_ivr(:,6),'--b')
xlim([0,t_final])
xline(S(3).TransientTime,'--r')
legend({'${\theta}$','$\hat{\theta}$'},'Interpreter','latex','Location','best')
xlabel('Time (s)')
ylabel('rad')
grid on

subplot(2,2,2)
plot(t_ivr,y_ivr(:,7),'r')
yline(double(u1_max),'--b')
yline(double(u1_min),'--b')
xlim([0,t_final])
legend('u_1')
ylabel('N')
grid on

subplot(2,2,4)
plot(t_ivr,y_ivr(:,8),'r')
yline(double(u2_max),'--b')
yline(double(u2_min),'--b')
xlim([0,t_final])
legend('u_2')
ylabel('N')
xlabel('Time (s)')
grid on

sgtitle('Initial Response for Outputs (left) & Inputs (right)');
set(f,"Position", [0 0 800 800])

In addition to the normal five subplot figure, I have also included the animation plot in the LiveScript as documentation of the performance. The als_f23_animate.p MATLAB p-code file can be found in the GitHub repository of the project.

als_f23_animate( t_ivr, y_ivr(:,1), y_ivr(:,2), y_ivr(:,3));

Animated GIF

For cases 6, 7, and 8, we will design and simulate another Output Feedback Controller (i.e. a state estimator with full state feedback control), but now we must also include Integral Control for reference tracking as defined in the notes.

We must expect to iterate on CASES 6, 7, and 8 until we achieve one acceptable design.

CASE 6. Design an LQI compensator

We must use LQR optimal control to design the state-feedback gains for this simulation. The MathWorks has provided a separate tool called LQI for computing the state-feedback gains when we are implementing state feedback with integral control. Note that LQI is identical to LQR except that LQI already knows that we are designing for reference tracking with an integral controller. We can NOT use LQR for this design, and we must use LQI! We must also use the LTI object for the open-loop system as input to the LQI function. Read the online help for LQI before starting this CASE!

First define the LQI cost function weighting matrices \(\mathbf{Q}\) and \(\mathbf{R}\). Note that LQI also accepts a weighting matrix \(\mathbf{N}\), which allows cross-weighting between states and controls, but we should NOT use this weighting matrix for this CASE! Make sure we understand the dimensions of each of these weighting matrices.

Q = diag([1 1 1 1 1 1 1 1 1])
Q = 9x9    
     1     0     0     0     0     0     0     0     0
     0     1     0     0     0     0     0     0     0
     0     0     1     0     0     0     0     0     0
     0     0     0     1     0     0     0     0     0
     0     0     0     0     1     0     0     0     0
     0     0     0     0     0     1     0     0     0
     0     0     0     0     0     0     1     0     0
     0     0     0     0     0     0     0     1     0
     0     0     0     0     0     0     0     0     1
R = diag([1 2])
R = 2x2    
     1     0
     0     2

Generate and display the optimal control gains and separate them into \({\mathbf{G}}_0\) and \({\mathbf{G}}_I\) matrices.

OL_int_A = [Aol, zeros(6,3);
    -Col, zeros(3,3)];
OL_int_B = [Bol;
    -Dol];
OL_int_sys = ss(OL_int_A,OL_int_B,[],[]);
[G,S,e]=lqi(OL_int_sys,Q,R)
G = 2x9    
   -3.2779   -2.4232    3.4413    1.1089   29.9386   11.5472    0.8981   -0.4397    0.0000
    1.1514    0.8555    5.0041    1.6054  -11.7328   -4.7892   -0.3109   -0.6351    0.0000

S = 9x9    
1.0e+08 *

    0.0507    0.0507    0.0000    0.0000   -0.0000   -0.0000   -0.0000   -0.0000   -0.4976
    0.0507    0.0507    0.0000    0.0000   -0.0000   -0.0000   -0.0000   -0.0000   -0.4976
    0.0000    0.0000    0.0000    0.0000   -0.0000   -0.0000   -0.0000   -0.0000    0.0000
    0.0000    0.0000    0.0000    0.0000   -0.0000   -0.0000   -0.0000   -0.0000    0.0000
   -0.0000   -0.0000   -0.0000   -0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
   -0.0000   -0.0000   -0.0000   -0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
   -0.0000   -0.0000   -0.0000   -0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
   -0.0000   -0.0000   -0.0000   -0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
   -0.4976   -0.4976    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    4.8816

e = 9x1 complex    
  -0.0000 + 0.0000i
  -0.1348 + 0.1327i
  -0.1348 - 0.1327i
  -0.5013 + 0.7953i
  -0.5013 - 0.7953i
  -0.9185 + 0.3250i
  -0.9185 - 0.3250i
  -1.6934 + 0.0000i
  -3.0023 + 0.0000i

CASE 7. Construct the OFC with Reference Inputs

Generate an augmented state-space system for the closed-loop Output Feedback Controller (state feedback with estimated states) and Reference Tracking. Our complete closed-loop state-space model must have:

{ State vector consisting of: \(\mathbf{x}\), \(\hat{\mathbf{x}}\), and \({\mathbf{x}}_I\) }

{ Output vector consisting of: \(\mathbf{y}\), \(\hat{\mathit{\mathbf{y}}}\), and \(\mathbf{u}\) }

{ Input vector consisting of: \(\mathbf{r}\) }

{ State Feedback Control law: \(u=-\left\lbrack \begin{array}{cc} {\mathbf{G}}_0 & {\mathbf{G}}_I \end{array}\right\rbrack \left\lbrack \begin{array}{c} \hat{\mathbf{x}} \\ {\mathbf{x}}_I \end{array}\right\rbrack\) }

After generating the closed-loop output feedback control state-space model, compute and display the closed-loop eigenvalues.

CASE 8. Simulate the Output Feedback Controller with Reference Inputs

For this CASE, we will simulate the closed-loop response of our control system using LSIM(); however, we must first generate a suitable reference profile as an input. Before we can generate a reference profile, we need to understand the initial position and the final desired position of the drone. These are new initial conditions, and a very different final position from any closed-loop simulation we have done before.

Initial Plant Conditions: The plant (drone) must start with the following initial states:

\(x\left(0\right)=-2\) m, \(y\left(0\right)=1\) m, and \(\theta \left(0\right)=0\) rad.

\(\dot{x} \left(0\right)=0\) m/s, \(\dot{y} \left(0\right)=0\) m/s, and \(\dot{\theta} \left(0\right)=0\) rad/s.

Initial Observer Conditions: We may assume that the initial observer states are identical to the initial plant states.

Final Plant Conditions: The plant (drone) must reach the following final states:

\(x\left(0\right)=+2\) m, \(y\left(0\right)=0\) m, and \(\theta \left(0\right)=0\) rad.

\(\dot{x} \left(0\right)=0\) m/s, \(\dot{y} \left(0\right)=0\) m/s, and \(\dot{\theta} \left(0\right)=0\) rad/s.

Evaluation Criteria:

  • Notice that the final state is on the opposite side of the 1.5m high red wall in the animation! The drone is NOT allowed to make contact with the wall at any point during the simulation.

  • The drone is NOT allowed to make contact with the ground at any point during the simulation except for the final time where the drone must sit on the ground.

    { The drone must reach the final state by 30 seconds (\(\textrm{tfinal}\le 30s\)). }

  • The control signals \(u_1\) and \(u_2\) must remain within the saturation bounds established in the Midterm Project.

    { There is NO formal settling time requirement; however, our actual output responses \(\mathbf{y}\) must track the reference responses \(\mathbf{r}\)reasonably accurately. For example, if our reference profiles have dwell regions, then there should be zero steady-state error after transients have subsided. If our reference profiles have ramp regions, then we should expect to see a steady-state error during the ramp, which looks like a delayed response. How closely we track the references will depend on our choice of weighting coefficients in the Q matrix. Note that a steady-state error during the ramp is perfectly acceptable as long as we choose the reference profile such that the drone does not make contact with the wall! }

  • For this simulation, the observer and the controller must both be enabled at t=0 (i.e. do NOT attempt a two-stage simulation like CASE 5).

  • We MUST use LSIM to simulate the closed-loop time response

    { Before constructing the reference profiles, we must first construct a time vector for the simulation. I recommend choosing a uniform step size of \(\textrm{ts}=0\ldotp 002\) seconds: }

\[ t=\textrm{ts}*\left\lbrack 0:\textrm{round}\left(\frac{\textrm{tfinal}}{\left.\textrm{ts}\right)}\right)\right\rbrack \]

  • For each time step in the time vector, we must generate a reference vector for each of the outputs. We will need to generate a particular path vs. time profile for the desired x and y position; however, I recommend choosing the following reference profile for the angular position:

\[ \theta_{\textrm{ref}} =\textrm{zeros}\left(\textrm{size}\left(t\right)\right) \]

  • I have used the same plot structure as before except that new signals that must be added. Each of the three output subplots on the left include the reference profile, the actual output, and the estimated output on the same axis.
  • Our main goal should be able to reach the final state in a minimum time.

SIMULATION RESULT

Animated GIF

Back to top
Source Code
---
title: "Full State Feedback Control Design"

description: "Optimal LQR control for drone landing using MATLAB, emphasizing strategic pole placement for precise descent"
categories: [Control]
image: Simulation_animation.gif
format:
    html: 
        code-fold: show
        code-overflow: wrap
        code-tools: true
---

# [ME5554/AOE5754/ECE5754 Applied Linear Systems - Project II]{style="color:rgb(213,80,0)"}

```matlab
clear;
set(0,'defaultlinelinewidth',2);
format short;
sympref('FloatingPointOutput',true);
```

# CASE 1. Prepare the Open-Loop Plant

The state-space LTI object we developed in the Project I for the open-loop plant should have been saved in a .mat file. 

<br>

Note: Before we proceed with the implementation, one should be responsible for correcting any errors they have in their open-loop plant! In this case, the plant is provided in my GitHub. One should NOT regenerate the open-loop LTI state-space plant in this LiveScript!!! The open-loop plant MUST be loaded from a .mat file (created in ME5554/AOE5754/ECE5754 Applied Linear Systems - Project I - Linearized LTI Modeling Feedback Control) to avoid discrepancies.

```matlab
load("MP_nishant_ss_model.mat")
```

After loading the .mat file, display the open-loop numerical state-space model. This model MUST display appropriately chosen statenames, inputnames (ok to use defaults here), and outputnames.

```matlab
OL_sys
```

``` text:output
OL_sys =
 
  A = 
                   X Position    X Velocity    Y Position    Y Velocity  Angular Posi  Angular Velo
   X Position               0             1             0             0             0             0
   X Velocity               0            -1             0             0         -9.81             0
   Y Position               0             0             0             1             0             0
   Y Velocity               0             0             0            -3             0             0
   Angular Posi             0             0             0             0             0             1
   Angular Velo             0             0             0             0             0          -1.6
 
  B = 
                 Varitional F  Varitional F
   X Position               0             0
   X Velocity               0             0
   Y Position               0             0
   Y Velocity             0.1           0.1
   Angular Posi             0             0
   Angular Velo        0.1302       -0.0898
 
  C = 
                   X Position    X Velocity    Y Position    Y Velocity  Angular Posi  Angular Velo
   X Position               1             0             0             0             0             0
   Y Position               0             0             1             0             0             0
   Angular Posi             0             0             0             0             1             0
 
  D = 
                 Varitional F  Varitional F
   X Position               0             0
   Y Position               0             0
   Angular Posi             0             0
 
Continuous-time state-space model.
Model Properties
```

Next, extract the Aol, Bol, Col, and Dol matrices, but do not display these results (because we already displayed them in the previous step).

```matlab
Aol = OL_sys.A;
Bol = OL_sys.B;
Col = OL_sys.C;
Dol = OL_sys.D;
```

Compute the open-loop eigenvalues (ol_poles), and the system dimensions (N, M, P). Assign these quantities to variables for use later in the LiveScript, but do not display them here.

```matlab
ol_poles = eig(Aol)
```

``` text:output
ol_poles = 6x1    
         0
   -1.0000
         0
   -3.0000
         0
   -1.6000
```

```matlab
N = size(Aol,2)
```

``` text:output
N = 6
```

```matlab
M = size(Bol,2)
```

``` text:output
M = 2
```

```matlab
P = size(Col,1)
```

``` text:output
P = 3
```

# CASE 2. Simulate Closed-Loop Control

We are given the following set of desired closed-loop poles which are known to result in 5% settling times (transient times) that are less than 3 seconds without saturating any of the control signals.

```matlab
cl_poles = [
    -29 + 15i
    -29 - 15i
    -32 + 0i
    -1.6 + 0.76i
    -1.6 - 0.76i
    -1.2 + 0i];
```

Use the `PLACE` function to compute the full state feedback control gain matrix, and display our result.

```matlab
G = place(Aol,Bol,cl_poles)
```

``` text:output
G = 2x6    
1.0e+03 *

   -1.1332   -0.8551    0.0362    0.1151    4.2932    0.2707
    2.4067    1.7537    0.4483    0.1343   -8.6992   -0.3186
```

Assume that the actual plant starts with the following initial conditions (NOTE: These are NEW IC's!!!):

```matlab
x0 = -2;            % meters
x0dot = 0.5;        % meters/sec
y0 = 1;             % meters
y0dot = 0;          % meters/sec
theta0 = -0.5;      % radians
theta0dot = 0.1;    % radians/sec
```

Assume that we have access to the full state vector, then construct a closed-loop state-space LTI object using the feedback gains from above.

```matlab
Acl = Aol - Bol*G
```

``` text:output
Acl = 6x6    
1.0e+03 *

         0    0.0010         0         0         0         0
         0   -0.0010         0         0   -0.0098         0
         0         0         0    0.0010         0         0
   -0.1273   -0.0899   -0.0484   -0.0279    0.4406    0.0048
         0         0         0         0         0    0.0010
    0.3637    0.2688    0.0355   -0.0029   -1.3402   -0.0655
```

```matlab
Ccl = Col - Dol*G
```

``` text:output
Ccl = 3x6    
     1     0     0     0     0     0
     0     0     1     0     0     0
     0     0     0     0     1     0
```

We must augment the output equation to output the original outputs as well as the control signals.

```matlab
Ccl = [Ccl; -G]
```

``` text:output
Ccl = 5x6    
1.0e+03 *

    0.0010         0         0         0         0         0
         0         0    0.0010         0         0         0
         0         0         0         0    0.0010         0
    1.1332    0.8551   -0.0362   -0.1151   -4.2932   -0.2707
   -2.4067   -1.7537   -0.4483   -0.1343    8.6992    0.3186
```

```matlab
CL_sys = ss(Acl,[],Ccl,[]);
```

Simulate the initial value response using `INITIAL`. Choose a final simulation time of 4 seconds for this simulation.

```matlab
t_final = 4; % seconds
xinit = [x0;x0dot;y0;y0dot;theta0;theta0dot];
[y_ivr,t_ivr] = initial(CL_sys,xinit,t_final);
```

Estimate the 5% settling times for each output.

```matlab
S = lsiminfo(y_ivr,t_ivr,"SettlingTimeThreshold",0.05);
S(1:3).TransientTime % Settling times in seconds
```

``` text:output
ans = 2.4134
ans = 2.6797
ans = 0.9707
```

Generate a plot like what was required for the Midterm Project with the three outputs on the left-hand side (separate axes for each output), and the control signals on the right-hand side of the figure (separate axes for each control signal). Annotate each plot with some kind of marker for the 5% settling times, the min and max allowable control signals, grid lines, axis labels, and units.

```matlab
f = figure();

u1_min = -2039;
u1_max = 1959;
u2_min = -2057.1;
u2_max = 1940.9;
outputNames = {'x', 'y', '\theta', 'u1','u2'};
outputUnits = {'m', 'm', 'rad','N', 'N'};
inputNames = {'u1', 'u2'};
inputUnits = {'N', 'N'};

for i = 1:3
    subplot(3, 2, -1+2*i); % odd numbers for left side
  
    plot(t_ivr, y_ivr(:, i));
    
    index = find((t_ivr-S(i).TransientTime)>0,1);
    if index<=size(t_ivr,1)
        xline(t_ivr(index),'--r')
    end
    clear index

    title(sprintf('%s', outputNames{i}));
    xlabel('Time (s)');
    ylabel(sprintf('%s', outputUnits{i}));
    grid on
end

for i = 1:2
    subplot(2, 2, i*2); % even numbers for right side
    
    plot(t_ivr, y_ivr(:, i+3)); % 4 and 5 for u1 and u2
    ylim([-2500,2500])
    if i == 1
        yline(double(u1_max),'--b')
        yline(double(u1_min),'--b')
    end
    if i == 2
        yline(double(u2_max),'--b')
        yline(double(u2_min),'--b')
    end

    title(sprintf('%s', inputNames{(i)}));
    xlabel('Time (s)');
    ylabel(sprintf('%s', inputUnits{(i)}));
    grid on
    
end
% Adjust the overall plot title
sgtitle('Initial Response for Outputs (left) & Inputs (right)');
set(f,"Position", [0 0 800 800])
```

![](Full_State_Feedback_Control_Design_media/figure_0.png)

In addition to the normal five subplot figure, we must also include the animation plot in the LiveScript as documentation of the performance.

```matlab
als_f23_animate( t_ivr, y_ivr(:,1), y_ivr(:,2), y_ivr(:,3));
```

# CASE 3. Evaluate Observability

First show that this system is completely observable using all three outputs.

```matlab
Ob = obsv(Aol,Col);
rank(Ob)
```

``` text:output
ans = 6
```

```matlab
N
```

``` text:output
N = 6
```

Rank of observability matrix is equal to number of states.

Then we determine whether this system is completely observable using any subset of outputs and justify our answers.

For individual outputs

```matlab
rank(obsv(Aol,Col(1,:)))
```

``` text:output
ans = 4
```

```matlab
rank(obsv(Aol,Col(2,:)))
```

``` text:output
ans = 2
```

```matlab
rank(obsv(Aol,Col(3,:)))
```

``` text:output
ans = 2
```

None of the individual outputs obs matrices are full rank

```matlab
rank(obsv(Aol,Col([1 2],:)))
```

``` text:output
ans = 6
```

```matlab
rank(obsv(Aol,Col([1 3],:)))
```

``` text:output
ans = 4
```

```matlab
rank(obsv(Aol,Col([2 3],:)))
```

``` text:output
ans = 4
```

We can see that the combination of ouput 1 & 2 (x & y), has a obs matrix which is full rank meaning the system is obervable using just there 2 outputs.

# CASE 4. Design a Full-Order Luenberger Observer

In this CASE, we will compute the gains for a state estimator using all three outputs. We must first choose appropriate pole locations for the state estimator (observer), and then compute and display the associated feedback gain matrix using the `PLACE` function.

```matlab
alpha = 1.5;
obs_poles = alpha*cl_poles;% alpha > 1 is faster observer than controller
K = place(Aol.',Col.',obs_poles).'
```

``` text:output
K = 6x3    
   47.0621   11.3921  -12.6548
   59.1867  -53.1150   13.2504
   -6.5959   78.7139  -41.1033
   63.5968    2.2363   -0.7294
   -5.7511   40.1340   10.2240
    1.9765   -0.8166    6.4001
```

To document our design, we must generate a plot that includes the open-loop poles from CASE 1, the closed-loop poles from CASE 2, and observer poles from this CASE. Use common markers for each set of poles (i.e. all open-loop poles must use the same marker, etc.) so that it is easy to distinguish the poles on the figure. Include a legend that only has three entries in it (open-loop, closed-loop, observer), corresponding to the three sets of poles. We must use the `SGRID` function (read the Matlab help!!) to generate appropriate grid lines. We will also annotate our plot with axis labels. Making sure our legend does not cover up any of the poles!

```matlab
f = figure;

wmax = round(max(abs([cl_poles;ol_poles;obs_poles])));
sgrid(0:0.1:0.9,0:wmax/5:wmax*1.5,"new")
h1 = plot(real(ol_poles),imag(ol_poles),'bx','markersize',10,'linewidth',2);
hold on
h2 = plot(real(cl_poles),imag(cl_poles),'gs','markersize',10,'linewidth',2);
h3 = plot(real(obs_poles),imag(obs_poles),'rd','markersize',10,'linewidth',2);
xlabel('Real (s)')
ylabel('Imag (s)')
xline(0,'m','linewidth',2)
yline(0,'m','linewidth',2)
legend([h1,h2,h3],'open loop','closed loop','estimator');
set(f,"Position", [0 0 600 600])
set(gca,'yaxislocation','right')
axis equal
```

![](Full_State_Feedback_Control_Design_media/figure_1.png)

# CASE 5. Simulate Output Feedback Control

For this CASE, we must simulate the complete output feedback control system, which will include the state estimator and full state feedback using the estimated states.

Using the equations for obersever followed by controller we can find, (L14/S29)

$$
\left\lbrack \begin{array}{c}
\dot{x} \\
\dot{\hat{x} } 
\end{array}\right\rbrack =\left\lbrack \begin{array}{cc}
A & -B\;G\\
K\;C & \hat{A} -K\hat{C} -\hat{B\;} G-K\left(D-\hat{D} \right)G
\end{array}\right\rbrack \left\lbrack \begin{array}{c}
x\\
\hat{x} 
\end{array}\right\rbrack
$$

$$
\left\lbrack \begin{array}{c}
y\\
\hat{y} 
\end{array}\right\rbrack =\left\lbrack \begin{array}{cc}
C & -D\;G\\
0 & \hat{C} -\hat{D} \;G
\end{array}\right\rbrack \left\lbrack \begin{array}{c}
x\\
\hat{x} 
\end{array}\right\rbrack
$$

Assumuing,

$$
\hat{A} =A\;\;\;\;\;\;\;\;\;\;\;\;\hat{B} =B\;\;\;\;\;\;\;\;\;\;\;\;\hat{C} =C\;\;\;\;\;\;\;\;\;\;\;\;\hat{D} =D
$$

We can simplify as,

$$
\left\lbrack \begin{array}{c}
\dot{x} \\
\dot{\hat{x} } 
\end{array}\right\rbrack =\left\lbrack \begin{array}{cc}
A & -B\;G\\
K\;C & A-K\;C-B\;G
\end{array}\right\rbrack \left\lbrack \begin{array}{c}
x\\
\hat{x} 
\end{array}\right\rbrack
$$

$$
\left\lbrack \begin{array}{c}
y\\
\hat{y} 
\end{array}\right\rbrack =\left\lbrack \begin{array}{cc}
C & -D\;G\\
0 & C-D\;G
\end{array}\right\rbrack \left\lbrack \begin{array}{c}
x\\
\hat{x} 
\end{array}\right\rbrack
$$

Augmenting to get back u to plot,

$$
\begin{bmatrix}
y \\
\hat{y} \\
u
\end{bmatrix}
=
\begin{bmatrix}
C & -DG \\
0 & C-DG \\
0 & -G
\end{bmatrix}
\begin{bmatrix}
x \\
\hat{x}
\end{bmatrix}
$$

In order to insure that the estimated states $\hat{\mathbf{x}}$ (and outputs $\hat{\mathbf{y}}$) have converged before control is enabled, we will first need to generate an augmented LTI closed-loop output-feedback control model with the state feedback gain $\mathbf{G}$ set to zero (it is easier to just use $\left(0*\mathbf{G}\right)$ in the matrices). This first LTI object will be used for the first **200 milliseconds** of the response. We must use the same initial states on the plant as defined in CASE 2 above, but the initial observer states must be zero. Our closed-loop state-space output equation must also be augmented to output the following three vectors: $\mathbf{y}$, $\hat{\mathbf{y}}$, and $\mathbf{u}$.

```matlab
xinit_OFC = [xinit;zeros(size(xinit))];

Aofc_G0 = [Aol, -Bol*G*0;
    K*Col, Aol-K*Col-Bol*G*0];
Cofc_G0 = [Col, -Dol*G*0;
    zeros(size(Col)), Col-Dol*G*0;
    zeros(size(G)), -G*0;
    eye(6,6),zeros(6,6);
    zeros(6,6),eye(6,6)];
OFC_sys_G0 = ss(Aofc_G0,[],Cofc_G0,[])
```

``` text:output
OFC_sys_G0 =
 
  A = 
             x1       x2       x3       x4       x5       x6       x7       x8       x9      x10      x11      x12
   x1         0        1        0        0        0        0        0        0        0        0        0        0
   x2         0       -1        0        0    -9.81        0        0        0        0        0        0        0
   x3         0        0        0        1        0        0        0        0        0        0        0        0
   x4         0        0        0       -3        0        0        0        0        0        0        0        0
   x5         0        0        0        0        0        1        0        0        0        0        0        0
   x6         0        0        0        0        0     -1.6        0        0        0        0        0        0
   x7     47.06        0    11.39        0   -12.65        0   -47.06        1   -11.39        0    12.65        0
   x8     59.19        0   -53.12        0    13.25        0   -59.19       -1    53.12        0   -23.06        0
   x9    -6.596        0    78.71        0    -41.1        0    6.596        0   -78.71        1     41.1        0
   x10     63.6        0    2.236        0  -0.7294        0    -63.6        0   -2.236       -3   0.7294        0
   x11   -5.751        0    40.13        0    10.22        0    5.751        0   -40.13        0   -10.22        1
   x12    1.976        0  -0.8166        0      6.4        0   -1.976        0   0.8166        0     -6.4     -1.6
 
  B = 
     Empty matrix: 12-by-0
 
  C = 
         x1   x2   x3   x4   x5   x6   x7   x8   x9  x10  x11  x12
   y1     1    0    0    0    0    0    0    0    0    0    0    0
   y2     0    0    1    0    0    0    0    0    0    0    0    0
   y3     0    0    0    0    1    0    0    0    0    0    0    0
   y4     0    0    0    0    0    0    1    0    0    0    0    0
   y5     0    0    0    0    0    0    0    0    1    0    0    0
   y6     0    0    0    0    0    0    0    0    0    0    1    0
   y7     0    0    0    0    0    0    0    0    0    0    0    0
   y8     0    0    0    0    0    0    0    0    0    0    0    0
   y9     1    0    0    0    0    0    0    0    0    0    0    0
   y10    0    1    0    0    0    0    0    0    0    0    0    0
   y11    0    0    1    0    0    0    0    0    0    0    0    0
   y12    0    0    0    1    0    0    0    0    0    0    0    0
   y13    0    0    0    0    1    0    0    0    0    0    0    0
   y14    0    0    0    0    0    1    0    0    0    0    0    0
   y15    0    0    0    0    0    0    1    0    0    0    0    0
   y16    0    0    0    0    0    0    0    1    0    0    0    0
   y17    0    0    0    0    0    0    0    0    1    0    0    0
   y18    0    0    0    0    0    0    0    0    0    1    0    0
   y19    0    0    0    0    0    0    0    0    0    0    1    0
   y20    0    0    0    0    0    0    0    0    0    0    0    1
 
  D = 
     Empty matrix: 20-by-0
 
Continuous-time state-space model.
Model Properties
```

Next, we must generate a second augmented LTI closed-loop output-feedback control model that is identical to the first model except that the state-feedback control gains are not zero, i.e. the controller is enabled. This second LTI object will be used for the remainder of the simulation.

```matlab
Aofc = [Aol, -Bol*G;
    K*Col, Aol-K*Col-Bol*G];
Cofc = [Col, -Dol*G;
    zeros(size(Col)), Col-Dol*G;
    zeros(size(G)), -G];
OFC_sys = ss(Aofc,[],Cofc,[])
```

``` text:output
OFC_sys =
 
  A = 
             x1       x2       x3       x4       x5       x6       x7       x8       x9      x10      x11      x12
   x1         0        1        0        0        0        0        0        0        0        0        0        0
   x2         0       -1        0        0    -9.81        0        0        0        0        0        0        0
   x3         0        0        0        1        0        0        0        0        0        0        0        0
   x4         0        0        0       -3        0        0   -127.3   -89.85   -48.45   -24.94    440.6    4.792
   x5         0        0        0        0        0        1        0        0        0        0        0        0
   x6         0        0        0        0        0     -1.6    363.7    268.8    35.54   -2.918    -1340   -63.86
   x7     47.06        0    11.39        0   -12.65        0   -47.06        1   -11.39        0    12.65        0
   x8     59.19        0   -53.12        0    13.25        0   -59.19       -1    53.12        0   -23.06        0
   x9    -6.596        0    78.71        0    -41.1        0    6.596        0   -78.71        1     41.1        0
   x10     63.6        0    2.236        0  -0.7294        0   -190.9   -89.85   -50.68   -27.94    441.3    4.792
   x11   -5.751        0    40.13        0    10.22        0    5.751        0   -40.13        0   -10.22        1
   x12    1.976        0  -0.8166        0      6.4        0    361.7    268.8    36.36   -2.918    -1347   -65.46
 
  B = 
     Empty matrix: 12-by-0
 
  C = 
           x1      x2      x3      x4      x5      x6      x7      x8      x9     x10     x11     x12
   y1       1       0       0       0       0       0       0       0       0       0       0       0
   y2       0       0       1       0       0       0       0       0       0       0       0       0
   y3       0       0       0       0       1       0       0       0       0       0       0       0
   y4       0       0       0       0       0       0       1       0       0       0       0       0
   y5       0       0       0       0       0       0       0       0       1       0       0       0
   y6       0       0       0       0       0       0       0       0       0       0       1       0
   y7       0       0       0       0       0       0    1133   855.1   -36.2  -115.1   -4293  -270.7
   y8       0       0       0       0       0       0   -2407   -1754  -448.3  -134.3    8699   318.6
 
  D = 
     Empty matrix: 8-by-0
 
Continuous-time state-space model.
Model Properties
```

We are now ready to simulate the closed-loop response, but it must happen in two separate simulation stages. In the first stage, we will use the INITIAL function to simulate for 200 milliseconds. At the end of this first simulation stage, we must capture the final state vector and use it as the initial state vector for the second simulation stage. The total simulation duration should be 4 seconds like above.

```matlab
t_break = 200e-03; % seconds
t_final = 4; % seconds
dt = 0.01;
t_1 = (0:dt:t_break);
t_2 = (t_break+dt):dt:t_final;
[y_1,t_1] = initial(OFC_sys_G0,xinit_OFC,t_1);
xinit_OFC_break = y_1(size(y_1,1),9:20).';
[y_2,t_2] = initial(OFC_sys, xinit_OFC_break,t_2);
y_ivr = [y_1(:,1:8);y_2(:,1:8)];
t_ivr = [t_1;t_2];
S = lsiminfo(y_ivr,t_ivr,"SettlingTimeThreshold",0.05);
S(1:3).TransientTime
```

``` text:output
ans = 2.5721
ans = 1.6099
ans = 2.1822
```

Generate the standard annotated plot (three output subplots on the left and two control signal subplots on the right), with appropriate annotations; but we must also plot the estimated output signals on top of the associated actual output signals. It must be clear from our plots that the estimated outputs all start from zero and converge to the actual outputs within 200 milliseconds. One must NOT try to get the observer states to converge too fast because it will be difficult to see the convergence on the plots. If we have done this step correctly, we should see 5% settling times that are comparable to the original design from CASE 2. The 5% settling times on the outputs must all be below 3 seconds, and we must not saturate either of the control signals.

```matlab
f = figure();
subplot(3,2,1)
plot(t_ivr,y_ivr(:,1),'r')
hold on
plot(t_ivr,y_ivr(:,4),'--b')
xline(S(1).TransientTime,'--r')
xlim([0,t_final])
legend({'${x}$','$\hat{x}$'},'Interpreter','latex','Location','best')
ylabel('m')
grid on

subplot(3,2,3)
plot(t_ivr,y_ivr(:,2),'g')
hold on
plot(t_ivr,y_ivr(:,5),'--black')
xlim([0,t_final])
xline(S(2).TransientTime,'--r')
legend({'${y}$','$\hat{y}$'},'Interpreter','latex','Location','best')
ylabel('m')
grid on

subplot(3,2,5)
plot(t_ivr,y_ivr(:,3),'m')
hold on
plot(t_ivr,y_ivr(:,6),'--b')
xlim([0,t_final])
xline(S(3).TransientTime,'--r')
legend({'${\theta}$','$\hat{\theta}$'},'Interpreter','latex','Location','best')
xlabel('Time (s)')
ylabel('rad')
grid on

subplot(2,2,2)
plot(t_ivr,y_ivr(:,7),'r')
yline(double(u1_max),'--b')
yline(double(u1_min),'--b')
xlim([0,t_final])
legend('u_1')
ylabel('N')
grid on

subplot(2,2,4)
plot(t_ivr,y_ivr(:,8),'r')
yline(double(u2_max),'--b')
yline(double(u2_min),'--b')
xlim([0,t_final])
legend('u_2')
ylabel('N')
xlabel('Time (s)')
grid on

sgtitle('Initial Response for Outputs (left) & Inputs (right)');
set(f,"Position", [0 0 800 800])
```

![](Full_State_Feedback_Control_Design_media/figure_2.png)

In addition to the normal five subplot figure, I have also included the animation plot in the LiveScript as documentation of the performance. The `als_f23_animate.p` MATLAB p-code file can be found in the GitHub repository of the project.

```matlab
als_f23_animate( t_ivr, y_ivr(:,1), y_ivr(:,2), y_ivr(:,3));
```

<img src="Full_State_Feedback_Control_Design_media/Simulation_animation.gif" alt="Animated GIF" width="50%"/>

## For cases 6, 7, and 8, we will design and simulate another Output Feedback Controller (i.e. a state estimator with full state feedback control), but now we must also include Integral Control for reference tracking as defined in the notes.

\matlabheadingtwo{This closed-loop simulation will require we to redesign the state-feedback control gains $\mathbf{G}$ (we cannot use the gains from CASE 2!); however, we should plan to use the same state estimator gains $\mathbf{K}$we designed in CASE 4.  we are welcome to redesign the state estimater gains (poles), but this is not recommended because it will add unnecessary time to our project completion.}

## We must expect to iterate on CASES 6, 7, and 8 until we achieve one acceptable design.

# CASE 6. Design an LQI compensator

We must use LQR optimal control to design the state-feedback gains for this simulation. The MathWorks has provided a separate tool called `LQI` for computing the state-feedback gains when we are implementing state feedback with integral control. Note that `LQI` is identical to `LQR` except that `LQI` already knows that we are designing for reference tracking with an integral controller. We can NOT use `LQR` for this design, and we must use `LQI!` We must also use the LTI object for the open-loop system as input to the `LQI` function. Read the online help for `LQI` before starting this CASE!

First define the LQI cost function weighting matrices $\mathbf{Q}$ and $\mathbf{R}$. Note that `LQI` also accepts a weighting matrix $\mathbf{N}$, which allows cross-weighting between states and controls, but we should NOT use this weighting matrix for this CASE! Make sure we understand the dimensions of each of these weighting matrices.

```matlab
Q = diag([1 1 1 1 1 1 1 1 1])
```

``` text:output
Q = 9x9    
     1     0     0     0     0     0     0     0     0
     0     1     0     0     0     0     0     0     0
     0     0     1     0     0     0     0     0     0
     0     0     0     1     0     0     0     0     0
     0     0     0     0     1     0     0     0     0
     0     0     0     0     0     1     0     0     0
     0     0     0     0     0     0     1     0     0
     0     0     0     0     0     0     0     1     0
     0     0     0     0     0     0     0     0     1
```

```matlab
R = diag([1 2])
```

``` text:output
R = 2x2    
     1     0
     0     2
```

Generate and display the optimal control gains and separate them into ${\mathbf{G}}_0$ and ${\mathbf{G}}_I$ matrices.

```matlab
OL_int_A = [Aol, zeros(6,3);
    -Col, zeros(3,3)];
OL_int_B = [Bol;
    -Dol];
OL_int_sys = ss(OL_int_A,OL_int_B,[],[]);
[G,S,e]=lqi(OL_int_sys,Q,R)
```

``` text:output
G = 2x9    
   -3.2779   -2.4232    3.4413    1.1089   29.9386   11.5472    0.8981   -0.4397    0.0000
    1.1514    0.8555    5.0041    1.6054  -11.7328   -4.7892   -0.3109   -0.6351    0.0000

S = 9x9    
1.0e+08 *

    0.0507    0.0507    0.0000    0.0000   -0.0000   -0.0000   -0.0000   -0.0000   -0.4976
    0.0507    0.0507    0.0000    0.0000   -0.0000   -0.0000   -0.0000   -0.0000   -0.4976
    0.0000    0.0000    0.0000    0.0000   -0.0000   -0.0000   -0.0000   -0.0000    0.0000
    0.0000    0.0000    0.0000    0.0000   -0.0000   -0.0000   -0.0000   -0.0000    0.0000
   -0.0000   -0.0000   -0.0000   -0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
   -0.0000   -0.0000   -0.0000   -0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
   -0.0000   -0.0000   -0.0000   -0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
   -0.0000   -0.0000   -0.0000   -0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
   -0.4976   -0.4976    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    4.8816

e = 9x1 complex    
  -0.0000 + 0.0000i
  -0.1348 + 0.1327i
  -0.1348 - 0.1327i
  -0.5013 + 0.7953i
  -0.5013 - 0.7953i
  -0.9185 + 0.3250i
  -0.9185 - 0.3250i
  -1.6934 + 0.0000i
  -3.0023 + 0.0000i
```

# CASE 7. Construct the OFC with Reference Inputs

Generate an augmented state-space system for the closed-loop Output Feedback Controller (state feedback with estimated states) and Reference Tracking. Our complete closed-loop state-space model must have:

\item

{ State vector consisting of: $\mathbf{x}$, $\hat{\mathbf{x}}$, and ${\mathbf{x}}_I$ }

\item

{ Output vector consisting of: $\mathbf{y}$, $\hat{\mathit{\mathbf{y}}}$, and $\mathbf{u}$ }

\item

{ Input vector consisting of: $\mathbf{r}$ }

\item

{ State Feedback Control law: $u=-\left\lbrack \begin{array}{cc} {\mathbf{G}}_0 & {\mathbf{G}}_I \end{array}\right\rbrack \left\lbrack \begin{array}{c} \hat{\mathbf{x}} \\ {\mathbf{x}}_I \end{array}\right\rbrack$ }

After generating the closed-loop output feedback control state-space model, compute and display the closed-loop eigenvalues.

```matlab
```

# CASE 8. Simulate the Output Feedback Controller with Reference Inputs

For this CASE, we will simulate the closed-loop response of our control system using `LSIM()`; however, we must first generate a suitable reference profile as an input. Before we can generate a reference profile, we need to understand the initial position and the final desired position of the drone. These are new initial conditions, and a very different final position from any closed-loop simulation we have done before.

**Initial Plant Conditions:** The plant (drone) must start with the following initial states:
         

$x\left(0\right)=-2$ m, $y\left(0\right)=1$ m, and $\theta \left(0\right)=0$ rad.

$\dot{x} \left(0\right)=0$ m/s, $\dot{y} \left(0\right)=0$ m/s, and $\dot{\theta} \left(0\right)=0$ rad/s.


**Initial Observer Conditions:** We may assume that the initial observer states are identical to the initial plant states.

**Final Plant Conditions:** The plant (drone) must reach the following final states:

        
$x\left(0\right)=+2$ m, $y\left(0\right)=0$ m, and $\theta \left(0\right)=0$ rad.

$\dot{x} \left(0\right)=0$ m/s, $\dot{y} \left(0\right)=0$ m/s, and $\dot{\theta} \left(0\right)=0$ rad/s.


**Evaluation Criteria:**

-   Notice that the final state is on the opposite side of the 1.5m high red wall in the animation! The drone is NOT allowed to make contact with the wall at any point during the simulation.\

-   The drone is NOT allowed to make contact with the ground at any point during the simulation except for the final time where the drone must sit on the ground.

    \item

    { The drone must reach the final state by 30 seconds ($\textrm{tfinal}\le 30s$). }

-   The control signals $u_1$ and $u_2$ must remain within the saturation bounds established in the Midterm Project.

    \item

    { There is NO formal settling time requirement; however, our actual output responses $\mathbf{y}$ must track the reference responses $\mathbf{r}$reasonably accurately. For example, if our reference profiles have dwell regions, then there should be zero steady-state error after transients have subsided. If our reference profiles have ramp regions, then we should expect to see a steady-state error during the ramp, which looks like a delayed response. How closely we track the references will depend on our choice of weighting coefficients in the Q matrix. Note that a steady-state error during the ramp is perfectly acceptable as long as we choose the reference profile such that the drone does not make contact with the wall! }

-   For this simulation, the observer and the controller must both be enabled at t=0 (i.e. do NOT attempt a two-stage simulation like CASE 5).

-   We MUST use `LSIM` to simulate the closed-loop time response

    \item

    { Before constructing the reference profiles, we must first construct a time vector for the simulation. I recommend choosing a uniform step size of $\textrm{ts}=0\ldotp 002$ seconds: }

$$
t=\textrm{ts}*\left\lbrack 0:\textrm{round}\left(\frac{\textrm{tfinal}}{\left.\textrm{ts}\right)}\right)\right\rbrack
$$

-   For each time step in the time vector, we must generate a reference vector for each of the outputs. We will need to generate a particular path vs. time profile for the desired x and y position; however, I recommend choosing the following reference profile for the angular position:

$$
\theta_{\textrm{ref}} =\textrm{zeros}\left(\textrm{size}\left(t\right)\right)
$$

-   I have used the same plot structure as before except that new signals that must be added. Each of the three output subplots on the left include the reference profile, the actual output, and the estimated output on the same axis.
-   Our main goal should be able to reach the final state in a minimum time.

## SIMULATION RESULT

<img src="Simulation_animation.gif" style="display: block; margin: 0 auto;" alt="Animated GIF" loop width="50%" height="100%">
Powered by Quarto
© 2024, Nishant Bharali

View source

Cookie Preferences
License: CC BY NC SA 4.0.