Main Content

このページは機械翻訳を使用して翻訳されました。元の英語を参照するには、ここをクリックします。

地上局アクセスの制約を満たしながら衛星群を最適化する

この例では、Global Optimization Toolbox ソルバーと Aerospace Toolbox ™ 関数を使用して、地上局のアクセス要件を満たす最小限の衛星コンステレーションを見つける方法を示します。この例では、4 つのステーションからなる地上局ネットワークが指定されていると想定しています。より複雑な最適化には、衛星だけでなく地上局の割り当ても含まれる可能性があります。アクセス要件は次のとおりです。

  • いくつかの関心地点は、少なくとも一定の割合の時間、衛星から観測可能でなければなりません。この例では、関心のあるポイントは地上局です。

  • 各衛星は、データをダウンロードするために、少なくとも一定の割合の時間、地上局を監視する必要があります。これらのパーセンテージは、accessPercentage (Aerospace Toolbox) 関数を使用して計算できます。

すべての衛星は同じ高度にあり、高度は最適化変数の 1 つです。衛星の軌道傾斜角も同じであり、これはもう 1 つの最適化変数です。衛星は複数の軌道面上にあり、最適化変数は 10 以下の整数です。最後に、軌道面あたりの衛星数は最適化変数であり、20 以下の整数です。

目的は衛星の総数を最小限に抑えることです。

星座の定義

ワークスペースに 4 つの地上局の座標を配置します。地上局を地図上にプロットします。

targets = table(...
    Size=[0,3],...
    VariableTypes=["string","double","double"],...
    VariableNames=["Name","Latitude","Longitude"]);

targets(1:4,:) = {...
    "Paris",        48.856614,      2.3522219  ;
    "Bruxelles",    50.850340,      4.351710   ;
    "Kiruna",       64.8557995,     20.2252821 ;    % Far from the equator
    "Kourou",       5.1597,         -52.6503  };    % Close to the equator
figure;
geoaxes("Basemap","satellite");
geoplot(targets.Latitude, targets.Longitude,"ko",MarkerFaceColor="red");
% Annotate the targets.
text(targets.Latitude,targets.Longitude,targets.Name,FontWeight="bold");

Figure contains an axes object with type geoaxes. The geoaxes object contains 5 objects of type line, text. One or more of the lines displays its values using only markers

地上局は、少なくとも 1 つの衛星によって、少なくとも minObservability = 90% の時間にわたって観測可能でなければなりません。この制約が星座に対して満たされているかどうかを計算するには、この例の最後にある observability ヘルパー関数を使用します。この関数は、衛星シナリオ関数 walkerDelta (Aerospace Toolbox) と Aerospace Toolbox で使用可能なアクセス メソッドを使用します。

minObservability = 0.9;

衛星シナリオを定義し、指定されたターゲットを地上局として含めます。シナリオとターゲット オブジェクトを mission という名前の構造体にバンドルします。

% Setup scenario
mission.Scenario = satelliteScenario(...
    datetime(2023,09,02,12,0,0),...
    datetime(2023,09,02,12,0,0) + days(1),...
    300);

% Ground stations
mission.Targets = groundStation(...
    mission.Scenario,...
    Name=targets.Name,...
    Latitude=targets.Latitude,...
    Longitude=targets.Longitude,...
    MinElevationAngle=30);

衛星コンステレーション定義の一部となるいくつかの量を定義します。詳細については、walkerDelta (Aerospace Toolbox) を参照してください。

phasing = 0;
argLat = 15; % deg

最適化問題の定義

この問題を最適化問題として定式化するには、最適化変数を定義し、次にこれらの変数に基づいて制約と目的関数を定義します。最適化問題を設定するための使いやすいインターフェースを提供する 問題ベースの最適化ワークフロー を使用します。

最適化変数

最適化問題には次の変数を使用します。

  • numPlanes — 軌道面の数、1 から 10 までの整数。

  • satPerPlane — 軌道面あたりの衛星数。1 から 20 までの整数。

  • inclination — 傾斜度数、40 から 80 までの実数スカラー。

  • altitude — メガメートル単位の高度、0.5 から 1.2 までの実数スカラー。この変数を他の変数とほぼ同じスケールにするには、一般的なキロメートルやメートルの代わりにメガメートルを使用します。この例の最後にある alt2radius ヘルパー関数は、高さをメガメートルからメートルに変換し、変換に地球の半径を含めます。

optimvar 関数を使用して最適化変数を定義します。変数の型 (整数または連続、記載されている場合) と境界を含めます。

numPlanes = optimvar("numPlanes",Type="integer",LowerBound=1,UpperBound=10);
satPerPlane = optimvar("satPerPlane",Type="integer",LowerBound=1,UpperBound=20);
inclination = optimvar("inclination",LowerBound=40,UpperBound=80);
altitude = optimvar("altitude",LowerBound=0.5,UpperBound=1.2);

目的関数

目的は、衛星の数を最小限に抑えることです。衛星の数とは、軌道面あたりの衛星の数と軌道面の数を掛け合わせたものです。

minimize  satPerPlane*numPlanes

この最小化の目的関数を含む最適化問題を作成します。

satelliteProb = optimproblem(Objective=satPerPlane*numPlanes,ObjectiveSense="minimize");

制約

この問題の制約は

visibilityminObservability

可視性は変数の複雑な関数です。この例の最後にある observability ヘルパー関数がこの量を計算します。このヘルパー関数を問題に配置するには、fcn2optimexpr 関数を使用して関数を最適化式に変換します。

visibility = fcn2optimexpr(@observability,numPlanes,satPerPlane,mission,...
                altitude,inclination,phasing,argLat);

最適化問題に制約を追加します。

satelliteProb.Constraints = visibility >= minObservability;

問題を確認します。

show(satelliteProb)
  OptimizationProblem : 

	Solve for:
       altitude, inclination, numPlanes, satPerPlane
	where:
       numPlanes, satPerPlane integer

	minimize :
       satPerPlane*numPlanes


	subject to :
       arg_LHS >= arg_RHS

       where:

             arg1 = observability(numPlanes, satPerPlane, extraParams{1}, altitude, inclination, 0, 15);
             arg_LHS = arg1;
             arg2 = 0.9;
             arg1 = arg2(ones(1,4));
             arg_RHS = arg1(:);

       extraParams

	variable bounds:
       0.5 <= altitude <= 1.2

       40 <= inclination <= 80

       1 <= numPlanes <= 10

       1 <= satPerPlane <= 20

問題を解く

問題の定式化は完了です。問題を解決するために使用できるソルバーを特定します。

[default,available] = solvers(satelliteProb)
default = 
"ga"
available = 1×3 string
    "ga"    "surrogateopt"    "gamultiobj"

このような時間のかかる整数制約の問題の場合、デフォルトの ga よりも surrogateopt ソルバーの方が効率的であることがよくあります。

surrogateopt ソルバーを使用して問題を解きます。ソリューションの進行状況を監視するには、optimplotfvalconstr プロット関数を指定します。より高速なソリューションを試すには、並列コンピューティングを指定します (Parallel Computing Toolbox™ が必要)。より優れた代替モデルを取得するには、MinSurrogatePoints オプションをデフォルト値の 2 倍の 40 に設定します。また、関数評価の最大数をデフォルト値より 100 多い 300 に設定します。並列で実行すると、surrogateopt は再現可能な結果を生成しません。したがって、実行環境がシリアルの場合は、再現可能な結果を得るためにランダム シードを設定してください。

rng default % For reproducibility
options = optimoptions("surrogateopt",PlotFcn=@optimplotfvalconstr,...
    UseParallel=true,MaxFunctionEvaluations=300,MinSurrogatePoints=40);
[sol,fval,exitflag,output] = solve(satelliteProb,Solver="surrogateopt",Options=options)
Solving problem using surrogateopt.

Figure Optimization Plot Function contains an axes object. The axes object with title Best Function Value: 112, xlabel Iteration, ylabel Function value contains 2 objects of type scatter. These objects represent Best function value, Best function value (infeasible).

surrogateopt stopped because it exceeded the function evaluation limit set by 
'options.MaxFunctionEvaluations'.
sol = struct with fields:
       altitude: 1.2000
    inclination: 62.8358
      numPlanes: 7
    satPerPlane: 16

fval = 112
exitflag = 
    SolverLimitExceeded

output = struct with fields:
        elapsedtime: 4.9507e+03
          funccount: 300
    constrviolation: 3.4602e-04
               ineq: [-0.1000 -0.1000 -0.1000 3.4602e-04]
           rngstate: [1×1 struct]
            message: 'surrogateopt stopped because it exceeded the function evaluation limit set by ↵'options.MaxFunctionEvaluations'.'
             solver: 'surrogateopt'

ディスプレイソリューション

最適化ソリューションから個別の変数を取得します。

numPlanes = sol.numPlanes
numPlanes = 7
numSatPerPlane = sol.satPerPlane
numSatPerPlane = 16
altitudeMm = sol.altitude
altitudeMm = 1.2000
inclination = sol.inclination
inclination = 62.8358
disp("Total number of satellites: " + evaluate(satelliteProb.Objective,sol))
Total number of satellites: 112

衛星群を視覚化します。

radiusMeters = alt2radius(altitudeMm);

mission.Constellation = walkerDelta(...
    mission.Scenario,...
    radiusMeters,...
    inclination,...
    numSatPerPlane*numPlanes,...
    numPlanes,...
    phasing,...
    ArgumentOfLatitude=argLat,...
    Name="Sat");

satelliteScenarioViewer(mission.Scenario);

constellation.png

補助関数

このコードは関数 alt2radius を作成します。

function radiusMeters = alt2radius(altitude)
% Convert altitude in Mm to radius in m.

% Add the radius of the Earth in Mm.
radius = altitude + 6.378137; 

% Convert the radius to meters.
radiusMeters = radius * 1000 * 1000;
end

このコードは関数 observability を作成します。

function targetVisibilityPct = observability(numPlanes,satPerPlane,mission,...
                altitude,inclination,phasing,argLat)
% The radius is in megameters (1000s of km) to keep the optimization search
% bounded in [0.5 1.5], which is roughly low Earth orbit.

% Convert the radius to meters for walkerDelta setup.
radiusMeters = alt2radius(altitude);

mission.Constellation = walkerDelta(...
    mission.Scenario,...
    radiusMeters,...
    inclination,...
    satPerPlane*numPlanes,...
    numPlanes,...
    phasing,...
    ArgumentOfLatitude=argLat,...
    Name="Sat");

% Compute the access for each target.
targetVisibilityPct = zeros(numel(mission.Targets),1);
for targetIdx = 1:numel(mission.Targets)
    accessAnalysis = access(mission.Constellation,mission.Targets(targetIdx));
    satAccess = accessStatus(accessAnalysis);
    systemAccess = any(satAccess,1);
    targetVisibilityPct(targetIdx) = mean(systemAccess); 
end
end

参考

関連するトピック