summaryrefslogtreecommitdiffstats
path: root/templates/lqr_maxPI.m
blob: 47ddd37a5143565da7334d5f2d98d0b17ff98fd3 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (c) 2023, Amon Lahr, Simon Muntwiler, Antoine Leeman & Fabian Flürenbrock Institute for Dynamic Systems and Control, ETH Zurich.
%
% All rights reserved.
%
% Please see the LICENSE file that has been included as part of this package.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% function [H, h] = lqr_maxPI(Q, R, params)
%     % Define the linear system
%     A = params.model.A;
%     B = params.model.B;
%     
%     % Define the state and input constraints
%     s_max=params.constraints.MaxAbsPositionXZ;
%     y_max=params.constraints.MaxAbsPositionY;
%     u_max = params.constraints.MaxAbsThrust;
%     
%     % Define the cost function
%     [K, ~, ~] = dlqr(A, B, Q, R);
%     P = dare(A, B, Q, R)
%     F = -inv((B'*P*B+R))*B'*P*A
%     A
%     B
%     Ax=[1,0;
%         0,1;
%         -1,0;
%         0,-1];
%     Ax=[1,0;
%         -1,0;
%         0,1;
%         0,-1];
%         Au=Ax;
%     bx=[s_max;s_max;s_max;s_max];
%     bu=[u_max;u_max;u_max;u_max];
%     H=[Ax;Au*F]/2
%     h=[bx;bu]
%     end

function [H, h] = lqr_maxPI(Q, R, params)
    % Define the linear system
    A = params.model.A;
    B = params.model.B;
    nx=params.model.nx;
    nu=params.model.nu;
    
    % Define the state and input constraints
    s_max= params.constraints.MaxAbsPositionXZ;
    y_max= params.constraints.MaxAbsPositionY;
    u_max = params.constraints.MaxAbsThrust;
    
    % Define the cost function
    K = -dlqr(A, B, Q, R);
    systemLQR = LTISystem('A', A+B*K);
%     absxmax = ones(nx,1)*s_max;
%     absumax = ones(nu,1)*u_max;
%     Xp = Polyhedron('A',[eye(nx); -eye(nx); K; -K], 'b', [absxmax;absxmax; absumax;absumax]);
    Hx=params.constraints.StateMatrix;
    hx=params.constraints.StateRHS;
    Hu=params.constraints.InputMatrix;
    hu=params.constraints.InputRHS;
    Xp = Polyhedron('A',[Hx;Hu*K], 'b', [hx;hu]);
    
    figure(1);
%     Xp.plot(),  alpha(0.25), title('Resulting State Constraints under LQR Control'), xlabel('x_1'), ylabel('x_2'), zlabel('x_3');
    systemLQR.x.with('setConstraint');
    systemLQR.x.setConstraint = Xp;
    InvSetLQR = systemLQR.invariantSet();
    InvSetLQR.plot(), alpha(0.25), title('Invariant Set for Triple Integrator under LQR Control'), xlabel('x_1'), ylabel('x_2'), zlabel('x_3');

    H=InvSetLQR.A
    h=InvSetLQR.b


end