Examine the "convergence rate" for an impulsive LP solution

As we increase the number of samples, the burn sequence converges to the exact solution. This happens slower at higher eccentricities, requiring a greater number of samples.

Details:
  Random initial and final states
  Vary the # of samples / orbit
  Vary eccentricity
  Vary duration from

Since version 7.
------------------------------------------------------------------------
Usage:
LPConvergenceDemo;
------------------------------------------------------------------------
See also AC, Mag, LPEccentric, FFEccGoals2Hills, OrbRate
------------------------------------------------------------------------

%-------------------------------------------------------------------------------
%   Copyright (c) 2004 Princeton Satellite Systems, Inc.
%   All rights reserved.
%-------------------------------------------------------------------------------

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% USER DEFINED DATA %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% semi-major axis
a    = 25000;

% eccentricity variation
e    = 0.1 : 0.1 : 0.3;

% maneuver duration
nOrb = 2;

% initial true anomaly
nu0  = 0;

% initial state
g0 = struct('y0',1,'xMax',1,'nu_xMax',0,'zMax',0,'nu_zMax',0);

% final state
gF = struct('y0',2,'xMax',1,'nu_xMax',0,'zMax',0,'nu_zMax',0);

% search data
changePercTol = 1;      % allowable percentage difference in solution
dNSPO         = 100;    % nSPO increment
f             = 1.5;
nSPO_e        = 1000;   % initial nSPO / e

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


n   = OrbRate(a);
nuF = nu0 + nOrb*2*pi;

nSFinal   = zeros(length(e),length(nOrb));
nSPOFinal = nSFinal;

for i=1:length(nOrb)

   for j=1:length(e)

      nSPO = round(nSPO_e*e(j));

      % compute initial/final states
      xH0 = FFEccGoals2Hills( e(j), nu0,    g0, n );
      xHF = FFEccGoals2Hills( e(j), nuF(i), gF, n );

      fprintf('\nnOrb = %f  ...  e = %f\n---------------------------\n',nOrb(i),e(j));
      oldDV = [];
      count = 0;

      % keep computing the control until it doesn't change much
      done = 0;
      while( ~done )

         count = count + 1;

         % update nS
         nS = round( nSPO*nOrb(i)^2 );

         % use LP to compute control sequence
         [aC,t,flag] = LPEccentric( e(j), n, xH0, xHF, nu0, nuF(i), nS );

         % compute the non-zero delta-v sequence
         k           = find(Mag(aC));
         dV          = Mag(aC) .* diff(t);
         newDV       = dV(k);

         % is the new delta-v sequence the same length as the previous one?
         if( length(oldDV) == length(newDV) )

            changePerc  = max( abs( (newDV-oldDV)/max(newDV) ) )*100;

            fprintf('%d - nSPO=%d, maximum percent change: %f\n',count,nSPO,changePerc);

            % if the percent change is within tolerance, we're done.
            % record nS and move on to next case.
            if( changePerc < changePercTol )
               done    = 1;
               nSFinal(i,j)   = nS;
               nSPOFinal(i,j) = nSPO;
            else
               %nSPO  = nSPO + dNSPO;
               nSPO  = round(nSPO*f);
               oldDV = newDV;
            end

         else

            fprintf('%d - nSPO=%d\n',count,nSPO);
            nSPO  = nSPO + dNSPO;
            oldDV = newDV;

         end

      end % while

   end % eccentricity loop

end % nOrb loop


%--------------------------------------

nOrb = 2.000000  ...  e = 0.100000
---------------------------
1 - nSPO=100
2 - nSPO=200, maximum percent change: 4.678900
3 - nSPO=300, maximum percent change: 2.144475
4 - nSPO=450, maximum percent change: 1.626262
5 - nSPO=675, maximum percent change: 0.706291

nOrb = 2.000000  ...  e = 0.200000
---------------------------
1 - nSPO=200
2 - nSPO=300, maximum percent change: 3.509740
3 - nSPO=450, maximum percent change: 1.700566
4 - nSPO=675, maximum percent change: 1.441010
5 - nSPO=1013, maximum percent change: 0.959928

nOrb = 2.000000  ...  e = 0.300000
---------------------------
1 - nSPO=300
2 - nSPO=400, maximum percent change: 2.890935
3 - nSPO=600, maximum percent change: 2.298576
4 - nSPO=900, maximum percent change: 1.781438
5 - nSPO=1350, maximum percent change: 1.307427
6 - nSPO=2025, maximum percent change: 0.769267

Published with MATLAB® R2019a