Relative orbit control using LQG
Note: This demo requies the Mu-Analysis and Synthesis Toolbox.
------------------------------------------------------------------------ See also AC, QCE, QCR, LinOrbLQG, DiscreteHills, LPCircular, LinOrb, OrbRate, Period ------------------------------------------------------------------------
Contents
------------------------------------------------------------------------- Copyright (c) 2004 Princeton Satellite Systems, Inc. All rights reserved. Since version 7. ------------------------------------------------------------------------- 24-Feb-2017: Update output from LinOrbLQG (4 matrices) -------------------------------------------------------------------------
if( ~exist('pck','file') | ~exist('sysic','file') ) msgbox('The mu-Analysis and Synthesis toolbox is required for this demo.','Missing software') return end % choose an altitude alt = 600; % mean orbit rate Re = 6378.14; sma = Re + alt; n = OrbRate( sma ); % open-loop plant [A,B] = LinOrb( [], n ); nI = size(B,2); nS = size(A,1); C = eye(nS); D = zeros(nS,nI); w.qE = eye(6); % State noise input matrix w.rE = eye(6); % Measurement noise covariance matrix w.q = (1e-4)*C'*C; % output weighting matrix w.r = (1e+4)*eye(3); % input weighting matrix KR = QCR( A, B, w.q, w.r ); % regulator (gain matrix) KE = QCE( A, C, w.qE, w.rE ); % estimator (gain matrix)
Combine the estimator and regulator into one controller
%-------------------------------------------------------- [aK,bK,cK,dK] = LinOrbLQG([sma,0,0,0,0,0],1e-4,1e4); % systems p = pck(A,B,C,D); % open-loop plant k = pck(aK,bK,cK,dK); % controller % interconnected system systemnames = ' p k '; inputvar = '[ r(6) ]'; outputvar = '[ p ; k ]'; input_to_p = '[ k ]'; input_to_k = '[ r ; p ]'; sysoutname = 'clp'; cleanupsysic = 'yes'; sysic; % initial & final state x0 = zeros(6,1); xF = [0;1;0;0;0;0]; % generate a trajectory to follow T = Period(sma); dTLP = 10; dur = round(T/dTLP)*dTLP; t1 = 0:dTLP:dur*2; aC = LPCircular(x0, xF, n, dur, dTLP ); aC = [aC, zeros(3,length(t1)-size(aC,2)-1)]; xS = DiscreteHills( x0, n, aC, dTLP ); r1 = xS; % choose a time-step for discretizing the system % and recompute the trajectory (linear interpolation btwn points) dT = 1; t2 = 0:dT:dur*2; r = interp1( t1, r1', t2 )'; % closed-loop system response tnew = 0:0.05:t2(end); rnew = interp1( t2, r', tnew )'; ref = vpck( rnew(:), tnew ); yCL = trsp(clp,ref,t2(end),.05,[xF;xF]*0); figure, vplot(sel(yCL,1:3,1)), title('position'); figure, vplot(sel(yCL,4:6,1)), title('velocity'); figure, vplot(sel(yCL,7:9,1)), title('control'); %-------------------------------------- % PSS internal file version information %--------------------------------------