Design and test the low frequency roll/yaw control system.

The script has three parts:

1. Investigate the orbit rate disturbance attenuation assuming that
   the orbit rate disturbance is caused by an inertially fixed torque.
2. Select gains and look at the torque transmission plots.
   The first set are closed loop.
   The second set are open loop (for reference purposes.)
3. Simulate the response to the orbit rate disturbance.
4. Must run ComStarDist.m first.
------------------------------------------------------------------------
See also Geo, ComStar, ComStarDist, RYDyn, SCHarm, SCSeries, C2DZOH,
MagPlot, TTPlots, Plot2D
------------------------------------------------------------------------------

Contents

%-------------------------------------------------------------------------------
%   Copyright (c) 1996 Princeton Satellite Systems, Inc.
%   All rights reserved.
%-------------------------------------------------------------------------------
%   Since version 2.
%-------------------------------------------------------------------------------

Constants

%----------
degToRad   = pi/180;
radToDeg   = 180/pi;

wo         = Geo;
hW         = ComStar('MWA Inertia')*ComStar('Nominal MWA Rate');
iAxis      = [1 3];
inr        = ComStar('MO Inertia');
hMWA       = ComStar('U MWA')*hW;
wO         = ComStar('Orbit Rate Unit Vector')*wo;
bMin       = 75e-9;

fprintf(1,'Inertia                        = %8.1f %8.1f %8.1f kg-m^2\n',inr(1,:));
fprintf(1,'                                 %8.1f %8.1f %8.1f kg-m^2\n',inr(2,:));
fprintf(1,'                                 %8.1f %8.1f %8.1f kg-m^2\n',inr(3,:));
fprintf(1,'MWA Momentum Vector            = %8.1f %8.1f %8.1f Nms\n',hMWA);
Inertia                        =   1640.1      0.0     -0.0 kg-m^2
                                      0.0    949.3      0.0 kg-m^2
                                     -0.0      0.0   2144.2 kg-m^2
MWA Momentum Vector            =      0.0    -62.8      0.0 Nms

Plant models

%-------------
[aP,bP,cP,dP,aL,bL,cL,dL] = RYDyn( inr, hMWA, wO, iAxis);

n    = 100;
nZ   = 6;
zeta = [0.1 0.2 0.5 0.7071 1.0 2.0];
k    = logspace(-2,1,n);
j    = sqrt(-1);
wN   = k*wo;
kX   = (wo^2-wN.^2)/wo;
f0   = 1.e-6/(wo*hW);
kZ   = zeros(nZ,n);
aX   = zeros(nZ,n);
aZ   = zeros(nZ,n);
for i = 1:nZ
  kZ(i,:) = 2*zeta(i)*wN;
  f       = f0./(j*kZ(i,:) - kX);
  aX(i,:) = abs( 2*wo.*f );
  aZ(i,:) = abs((j*wo + kZ(i,:) - j*(kX-wo)).*f);
end

Plot2D(k,aX*180/pi,'wN/wo','Roll (deg/micro Nm)','Roll Closed Loop Gain','xlog')
Plot2D(k,aZ*180/pi,'wN/wo','Yaw (deg/micro Nm)','Yaw Closed Loop Gain','xlog')
Plot2D(k,kZ,'wN/wo','kZ','Roll Angle to Yaw Torque Gain','xlog')
Plot2D(k,kX,'wN/wo','kX','Roll Angle to Roll Torque Gain','xlog')

Skew dipole control system simulation

%--------------------------------------
zeta  = 2.5;
wN    = 5*wo;
kRY   = [(wo^2-wN^2)/wo;2*zeta*wN]; % These gains attenuate the roll/yaw dc disturbance by 0.5
aCL   = aL + bL*kRY*[1 0]*hW;
dT    = 864;
s     = eig(aCL);
fprintf(1,'Closed loop eigenvalues        = %12.4e ± j%12.4e (rad/sec)\n',real(s(1)), imag(s(1)));
fprintf(1,'Roll angle to roll torque gain = %12.8f (Nm/deg)\n',kRY(1)*hW*pi/180);
fprintf(1,'Roll angle to yaw torque gain  = %12.8f (Nm/deg)\n',kRY(2)*hW*pi/180);

[a,b] = C2DZOH(aCL,bL,dT);
Closed loop eigenvalues        =  -1.7469e-03 ± j  0.0000e+00 (rad/sec)
Roll angle to roll torque gain =  -0.00191921 (Nm/deg)
Roll angle to yaw torque gain  =   0.00199917 (Nm/deg)

The torque transmission plots

%------------------------------
titlesRY = [
            'Roll Torque to Roll Angle';...
            'Roll Torque to Yaw Angle ';...
            'Yaw Torque to Roll Angle ';...
            'Yaw Torque to Yaw Angle  ';...
            ];
wP        = logspace(-6,-2,300);
mag       = MagPlot(aCL,bL,cL,dL,1:2,1:2,wP);
TTPlots(radToDeg*1.e-6*mag,wP,titlesRY,'deg/microNm','Closed Loop')
[mag,io] = MagPlot(aL,bL,cL,dL,1:2,1:2,wP);
TTPlots(radToDeg*1.e-6*mag,wP,titlesRY,'deg/microNm','Open Loop')

load DistModl

tHarm  = tHarm([1 3],:);

x      = [0;0];
t      = 0;
nSim   = 1000;
xPlot  = zeros(length(x),nSim);
uPlot  = zeros(2,nSim);
for k = 1:nSim
  t          = t + dT;
  [s, c]     = SCHarm( wo*t, 6 );
  u          = tHarm(:,1) + tHarm(:,2:7)*s' + tHarm(:,8:13)*c';
  x          = a*x + b*u*1.e-6;
  xPlot(:,k) = x;
  uPlot(:,k) = kRY*x(1)*hW/bMin;
end

rYGain = kRY*hW;

pData = fileparts(mfilename('fullpath'));
save(fullfile(pData,'RYC'),'rYGain', 'aP', 'bP', 'cP', 'dP' );

Plot2D( (0:(nSim-1))*dT,xPlot*180/pi,'Time (sec)','Roll/yaw (deg)','Solstice Roll/Yaw Errors')
Plot2D( (0:(nSim-1))*dT,uPlot,'Time (sec)','Dipole (ATM^2)','Solstice Torquer Demand')

Figui

maxATM = max(abs(uPlot'));

fprintf(1,'Maximum Roll Dipole = %12.8f (ATM^2)\n',maxATM(1));
fprintf(1,'Maximum Yaw Dipole  = %12.8f (ATM^2)\n',maxATM(2));
Maximum Roll Dipole = 405.21831149 (ATM^2)
Maximum Yaw Dipole  = 422.10240780 (ATM^2)

Harmonic expansion

%-------------------
t = (0:(nSim-1))*dT;

tHarm = [ SCSeries( t,6,6,wo,xPlot(1,:)*180/pi )';...
          SCSeries( t,6,6,wo,xPlot(2,:)*180/pi )'];

kPrint = [1:4 8:10];
fprintf(1,'       Bias    Sin(wo)   Sin(2*wo)  Sin(3*wo)    Cos(wo)  Cos(2*wo)  Cos(3*wo)\n');
fprintf(1,'x %9.4f  %9.4f  %9.4f  %9.4f  %9.4f  %9.4f  %9.4f     deg\n',  tHarm(1, kPrint) );
fprintf(1,'z %9.4f  %9.4f  %9.4f  %9.4f  %9.4f  %9.4f  %9.4f     deg\n\n',tHarm(2, kPrint) );


%--------------------------------------
       Bias    Sin(wo)   Sin(2*wo)  Sin(3*wo)    Cos(wo)  Cos(2*wo)  Cos(3*wo)
x    0.0030     0.0083     0.0052     0.0001     0.0006     0.0031    -0.0000     deg
z    0.0753     0.0904    -0.0218    -0.0017     0.1149     0.0683     0.0003     deg