Manually compute LQI() gains using Algebraic Riccati Equation / Hamiltonian
23 次查看(过去 30 天)
显示 更早的评论
I'm attempting to compile a Matlab program for embedded use, and that program also includes lqr() and lqi(). Since neither are compilable functions for compiling, I'm solving LQR and LQI manually, ultimately hoping to end up only with functions that can be compiled.
I'm doing this via two Hamiltonian methods, with pole placement or to find the Algebraic Riccati Equation solution P. This is working correctly (solutions P and gains K match lqr() outputs) for a continuous-time basic plant without integrator. But, for the discrete case it's failing for both methods I try, with results that don't match Matlab's results (i won't go into the discrete implementation in detail, for clarity, unless asked)
My question is:
- Is there a better manual method to finding lqi gains, especially one that generalizes to the discrete case more easily? I'm having trouble using the below methods for the discrete case, even with needed changes (see links in this post) to handle the discrete case. I can create a separate post for troubleshooting that, but for now, i'm just curious if there are simpler manual methods -- or built-in methods that Matlab allows for compiling into executables.
- If I set Q = eye(3) for LQI case, it's suspicious that all lqi gains are 1's, and that lqr() and lqi(1:2) gains don't match. Why would gains be almost exactly 1 from Q = eye(3)?
Details
I'm using two manual processes. One uses poles placement via acker(), and the other works via Schur decomposition via schur(), as seen in this link.
Both methods work for LQR (plant without integral action). But neither works with LQI (plant with integral action), insofar as matching the base lqr() gains. Of course, the plant with integral action won't use the integral action in open-loop, so it's more a structural addition.
Method 1: Pole placement method
Method 2: Schur method
Code
clc
close all
A = [0 1;
-4e5 -400];
B = [0; 4e5];
C = [1 0];
sys = ss(A, B, C, 0);
% i denotes "integral"
Ai = [A, zeros(2,1); -C, 0];
Bi = [B; 0];
Ci = [C 0];
sys_aug = ss(Ai, Bi, Ci, 0);
% ========================
% LQR
% x = [x1; x2]
Q = eye(size(A));
Q(1,1) = 1e5;
Q(2,2) = 0.1;
R = 1;
K_lqr = lqr(sys, Q, R)
% ~~~ Method 1 ~~~
hamiltonian = [A, -B/R*B';
-Q, -A'];
eigs_H = eig(hamiltonian);
whichNegEig = find(real(eigs_H)<0);
% This matches K_lqr
K_hamilt = acker(A, B, eigs_H(whichNegEig))
% ~~~ Method 2 ~~~
[U, S] = schur(hamiltonian);
[U, S] = ordschur(U, S, 'lhp');
% Find U21 and U11: U is a square matrix.
[m,n] = size(U);
U11 = U(1:(m/2), 1:(n/2));
U21 = U((m/2+1):m, 1:(n/2));
% Find solution P to A.R.E.
P = U21 * inv(U11);
% Double-check solution using built-in care()
P_care = care(A, B, Q, R) ;
% This matches K_lqr
K_schur = inv(R) * B' * P
Acl_lqr_manual = [sys.A - sys.B * K_schur];
syscl_lqr_manual = ss(Acl_lqr_manual, B* K_lqr * [1;0], C, 0);
figure;
step(syscl_lqr_manual)
hold on;
% ========================
% LQI: Integral action
% x = [x1; x2; xi]
Qi = eye(size(Ai));
Qi(1,1) = Q(1,1);
Qi(2,2) = Q(2,2);
Qi(3,3) = 1e12;
Ri = 1;
K_lqr_i = lqr(sys_aug, Qi, Ri)
K_lqi_i2 = lqi(sys, Qi, Ri) % just double-check
% ~~~ Method 1 ~~~
hamiltonian_i = [Ai, -Bi/Ri*Bi';
-Qi, -Ai'];
eigs_Hi = eig(hamiltonian_i);
whichNegEig_i = find(real(eigs_Hi)<0);
% This is all NaN due to 0 column
K_hamilt_i = acker(Ai, Bi, eigs_Hi(whichNegEig_i))
% ~~~ Method 2 ~~~
[Ui, Si] = schur(hamiltonian_i);
[Ui, Si] = ordschur(Ui, Si, 'lhp');
% Find U21 and U11: U is a square matrix.
[m,n] = size(Ui);
U11i = Ui(1:(m/2), 1:(n/2));
U21i = Ui((m/2+1):m, 1:(n/2));
% Find solution P to A.R.E.
Pi = U21i * inv(U11i);
% Double-check solution using built-in care()
Pi_care = care(Ai, Bi, Qi, Ri);
% This ~matches K_lqr_i, but not K_lqi_i.
K_schur_i = inv(Ri) * Bi' * Pi
Acl_lqi_manual = [sys.A - sys.B * K_schur_i(1:2), -sys.B * K_schur_i(3);
-sys.C, 0];
% Ref is through integ K(Ki*xi - x), ie
% Ref feedthrough only to xi, ie K(r-x), vs above w/ ref through integ K(Ki*xi - x)
Bcl_lqi_manual = [0; 0; 1];
Ccl_lqi_manual = [sys.C 0];
syscl_lqi_manual = ss(Acl_lqi_manual, Bcl_lqi_manual, Ccl_lqi_manual, 0);
step(syscl_lqi_manual)
step(sys)
legend('manual lqr', ' manual lqi', 'OL')
0 个评论
采纳的回答
Paul
2023-7-11
Hi John,
I just copy/pasted the code from the question (not the attached .m file) with some semicolons to suppress output and it seems to run fine. The gains for the plant all match when compute with the different methods and the gains for the augmented plant all match when computed with the different methods.
Perhaps I don't understand the concern with results from the code. Are you expecting some sort of commonality between the values of K_lqr and K_lqr_i?
Also, consider using place instead of acker, at least if the closed loop poles don't have multiplicity greater than one.
format long e
A = [0 1;
-4e5 -400];
B = [0; 4e5];
C = [1 0];
sys = ss(A, B, C, 0);
% i denotes "integral"
Ai = [A, zeros(2,1); -C, 0];
Bi = [B; 0];
Ci = [C 0];
sys_aug = ss(Ai, Bi, Ci, 0);
%% ========================
% LQR
% x = [x1; x2]
Q = eye(size(A));
R = 1;
K_lqr = lqr(sys, Q, R)
% ~~~ Method 1 ~~~
hamiltonian = [A, -B/R*B';
-Q, -A'];
eigs_H = eig(hamiltonian);
whichNegEig = find(real(eigs_H)<0);
% This matches K_lqr
K_hamilt = acker(A, B, eigs_H(whichNegEig))
% ~~~ Method 2 ~~~
[U, S] = schur(hamiltonian);
[U, S] = ordschur(U, S, 'lhp');
% Find U21 and U11: U is a square matrix.
[m,n] = size(U);
U11 = U(1:(m/2), 1:(n/2));
U21 = U((m/2+1):m, 1:(n/2));
% Find solution P to A.R.E.
P = U21 * inv(U11);
% Double-check solution using built-in care()
P_care = care(A, B, Q, R);
% This matches K_lqr
K_schur = inv(R) * B' * P
%% ========================
% LQI: Integral action
% x = [x1; x2; xi]
Qi = eye(size(Ai));
Ri = 1;
K_lqr_i = lqr(sys_aug, Qi, Ri)
K_lqi_i = lqi(sys, Qi, Ri) % just double-check
% ~~~ Method 1 ~~~
hamiltonian_i = [Ai, -Bi/Ri*Bi';
-Qi, -Ai'];
eigs_Hi = eig(hamiltonian_i);
whichNegEig_i = find(real(eigs_Hi)<0);
% This is all NaN due to 0 column
K_hamilt_i = acker(Ai, Bi, eigs_Hi(whichNegEig_i))
% ~~~ Method 2 ~~~
[Ui, Si] = schur(hamiltonian_i);
[Ui, Si] = ordschur(Ui, Si, 'lhp');
% Find U21 and U11: U is a square matrix.
[m,n] = size(Ui);
U11i = Ui(1:(m/2), 1:(n/2));
U21i = Ui((m/2+1):m, 1:(n/2));
% Find solution P to A.R.E.
Pi = U21i * inv(U11i);
% Double-check solution using built-in care()
Pi_care = care(Ai, Bi, Qi, Ri);
% This ~matches K_lqr_i, but not K_lqi_i.
K_schur_i = inv(Ri) * Bi' * Pi
5 个评论
Paul
2023-7-20
移动:Walter Roberson
2023-8-8
I've read this comment a few times and I don't fully understand the questions. I'll try to answer anyway basedon what I think the questions are poking at.
Suppose we have a continuous plant, modeled with Ac, Bc and we want to design an LQR controller that will be implemented in a discrete-time compensator.
Option 1: Design Kc_lqr using lqr based on Qc and Rc, where Qc and Rc define the continuous-time quadaratic cost integral. Usually this will involve some iteration on Qc and Rc until a satisfactory design is achieved. Because Kc_lqr is just a gain matrix, implement it direactly as a gain matrix in the discrete-time compensator implementation (just like one would do with any classical proportional control). Note that there is no guarantee that Kc_lqr is "optimal" in the mathematical sense when combined with the discrete-time approximation of the plant.
Option 2. Discretize the plant model using c2d (typically with zoh option) to get Ad and Bd. Design Kd_lqr using lqr or dlqr based on Qd and Rd, where Qd and Rd define the discrete-time quadratic cost sum. Usually this will involve some iteration on Qd and Rd until a satisfactory design is achieved. Implement Kd_lqr directly in discrete-time compensation. Kd_lqr is "optimal" for the discrete time approximation of the plant relative to the quadratic cost sum.
Option 3: Start with option 1, design K based on Ac, Bc, Qc and Rc. When the continuous-time design is satisfactory compute the discretized plant model Ad,Bd from Ac,Bc. Compute Qd,Rd from Qc,Rc so that the discrete-time quadratic cost sum in some sense approximates the continuous-time quadratic cost integral. This transformation from Qc,Rc to Qd,Rd can be done approximately or "exactly" (quotation marks because there is never anything exact going from continous to discrete); find a good text book or course notes to see how (or look at the Algorithms section of lqrd. Design Kd_lqr based on Ad, Bd, Qd, Rd using lqr or dlqr. Implement Kd_lqr in the discrete-time compensation. The idea here is to compute Kd_lqr that minimizes the discrete-time quadratic cost (computed from the discrete-time plant approximation), which in turn is close to "equivalent" in some sense as the continuous-time quadratic cost.
Option 4: Start with option 1, design K using lqr based on Ac, Bc, Qc and Rc. When the the continuous-time design is satisfactory compute Kd_lqrd using lqrd based on Ac, Bc, Qc and Rc. This option is doing the heavy lifting to compute the Qd, Rd (and Nd) that would be the ideal approach in Option 3.
Option 5: I suppose it's feasible to try iterate on Qc and Rc as inputs to lqrd to compute a Kd_lqrd until satisfactory design is achieved. I don't see any advantages to this approach.
更多回答(0 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!