このページの翻訳は最新ではありません。ここをクリックして、英語の最新版を参照してください。

離散化された最適軌跡 (問題ベース)

この例では、問題ベースのアプローチを使用して、離散化された最適軌跡問題を解く方法を示します。軌跡の重力は一定とし、印加される力には限度があり、空気抵抗は考慮しません。この解法プロセスでは、固定期間 $$T$$ の間の最適な軌跡を求め、その解を使用して最適な $$T$$、つまりコストを最小限にする時間を見つけます。最後のセクションでは、空気抵抗を導入する方法を示します。

最適化変数を使用して ODE のパラメーターを当てはめるの例のような離散化されていない最適化と比較して、離散化された最適化は常微分方程式 (ODE) を解く場合正確ではありません。しかし、離散化された最適化は、シミュレーションまたは常微分方程式の最適化の説明にあるとおり、可変ステップ ODE ソルバーに含まれるノイズの影響を受けません。また、離散化された最適化の方がカスタマイズしやすくなることがあり、モデル化も簡単です。最後に、離散化された最適化では、最適化プロセスで自動微分を活用できます。問題ベースの最適化における自動微分の効果を参照してください。

問題の説明

問題は、制御されたジェット機を使用して、時間 $$0$$ における位置 $$p_0$$ から時間 $$T$$ における位置 $$p_F$$ までオブジェクトを移動することです。位置はベクトル $$p(t)$$、速度はベクトル $$v(t)$$、印加された加速度はベクトル $$a(t)$$ で表します。連続時間では、重力を含めた運動方程式は次のようになります。

$$\frac{dp}{dt}=v(t)$$

$$\frac{dv}{dt}=a(t)+g*[0,0,-1]$$.

時間を離散化して問題を解きます。時間をサイズ $$t=T/N$$ で $$N$$ 個に等分割します。タイム ステップ $$i$$ における位置はベクトル $$p(i)$$、速度はベクトル $$v(i)$$、印加された加速度はベクトル $$a(i)$$ となります。ODE モデルをかなり正確に表現する、一連の方程式を作成できます。近似運動方程式は次のようになります。

$$v(i)=v(i-1)+t(a(i-1)+g)$$

$$\begin{array}{l}p(i)=p(i-1)+t\left(\frac{v(i-1)+v(i)}{2}\right)\\ =p(i-1)+tv(i-1)+t^2\frac{a(i-1)+g}{2}.\end{array}$$

前述の方程式では、2 点 (台形則) 近似を使用して速度を積分し、1 点 (オイラー) 近似を使用して加速度を積分します。この積分法で得られる方程式は簡単なもので、ステップ $$i$$ における位置および速度は、ステップ $$i-1$$ における位置、速度、および加速度のみに依存しています。この方程式は、空気抵抗を考慮するように変更するのも容易です。

境界条件は、$$p(1)=p_0$$、$$p(N)=p_F$$、および $$v(1)=v(N)=[0,0,0]$$ です。初期位置および最終位置は次のように設定します。

p0 = [0 0 0];
pF = [5 10 3];

ジェット機を使用して、時間 $$t$$ の間に力 $$a$$ を発生させるのにかかるコストは、$$\Vert a\Vert t$$ です。総コストは、加速度と $$t$$ の積のノルムの和です。

$$Cost=\sum _{i=1}^{N-1}\Vert a(i)\Vert t.$$

このコストを最適化変数の線形コストに変換するには、変数 $$s(i)$$ を作成し、関連する 2 次錐制約を作成します。

$$\begin{array}{l}Cost=\sum _{i=1}^{N-1}s(i)t\\ \Vert a(i)\Vert \le s(i).\end{array}$$

すべてのタイム ステップについて、加速度のノルムが定数 Amax によって制限される追加制約を課します。

$$\Vert a(i)\Vert \le Amax.$$

これらの制約も 2 次錐制約です。制約が線形または 2 次錐制約であり、目的関数が線形であるため、solve は coneprog ソルバーを呼び出して問題を解きます。

次のコードは、固定の時間 T の最適化問題を作成します。このコードでは、運動方程式を問題の制約として組み込んでいます。ライブ スクリプトを使用してこの例を実行することで、setupproblem.m 関数ファイルにアクセスできます。この関数には空気抵抗の引数が含まれています。空気抵抗を考慮するモデルの場合は、air = true と設定してください。空気抵抗の定義については、空気抵抗の導入のセクションを参照してください。

type setupproblem

function trajectoryprob = setupproblem(T,air)
    if nargin == 1
        air = false;
    end
    N = 50;
    g = [0 0 -9.81];
    p0 = [0 0 0];
    pF = [5 10 3];
    Amax = 25;
    t = T/N;
    p = optimvar("p",N,3);
    v = optimvar("v",N,3);
    a = optimvar("a",N-1,3,"LowerBound",-Amax,"UpperBound",Amax);
    trajectoryprob = optimproblem;
    s = optimvar("s",N-1,"LowerBound",0,"UpperBound",3*Amax);
    trajectoryprob.Objective = sum(s)*t;
    scons = optimconstr(N-1);
    for i = 1:(N-1)
        scons(i) = norm(a(i,:)) <= s(i);
    end
    acons = optimconstr(N-1);
    for i = 1:(N-1)
        acons(i) = norm(a(i,:)) <= Amax;
    end
    vcons = optimconstr(N+1,3);
    vcons(1,:) = v(1,:) == [0 0 0];
    if air
        vcons(2:N,:) = v(2:N,:) == v(1:(N-1),:)*exp(-t) + t*(a(1:(N-1),:) + repmat(g,N-1,1));
    else
        vcons(2:N,:) = v(2:N,:) == v(1:(N-1),:) + t*(a(1:(N-1),:) + repmat(g,N-1,1));
    end
    vcons(N+1,:) = v(N,:) == [0 0 0];
    pcons = optimconstr(N+1,3);
    pcons(1,:) = p(1,:) == p0;
    if air
        pcons(2:N,:) = p(2:N,:) == p(1:(N-1),:) + t*(1+exp(-t))/2*v(1:(N-1),:) + t^2/2*(a(1:(N-1),:) + repmat(g,N-1,1));
    else
        pcons(2:N,:) = p(2:N,:) == p(1:(N-1),:) + t*v(1:(N-1),:) + t^2/2*(a(1:(N-1),:) + repmat(g,N-1,1));
    end
    pcons((N+1),:) = p(N,:) == pF;
    trajectoryprob.Constraints.acons = acons;
    trajectoryprob.Constraints.scons = scons;
    trajectoryprob.Constraints.vcons = vcons;
    trajectoryprob.Constraints.pcons = pcons;
end

T = 20 について問題を解く

時間 $$T=20$$ の軌跡問題を作成し、解きます。

trajprob = setupproblem(20);
[sol,fval,eflag,output] = solve(trajprob)

Solving problem using coneprog.
Optimal solution found.

sol = struct with fields:
    a: [49x3 double]
    p: [50x3 double]
    s: [49x1 double]
    v: [50x3 double]

fval = 192.2989

eflag = 
    OptimalSolution

output = struct with fields:
           iterations: 8
    primalfeasibility: 3.2932e-07
      dualfeasibility: 2.9508e-07
           dualitygap: 1.7343e-08
            algorithm: 'interior-point'
         linearsolver: 'prodchol'
              message: 'Optimal solution found.'
               solver: 'coneprog'

この例の終わりに掲載されている補助関数 plottrajandaccel を呼び出して、軌跡および加速度のノルムをプロットします。

plottrajandaccel(sol,p0,pF)

加速度は、やや低下している終点付近を除いて、常に重力加速度 (9.81) 付近にあります。

最小コストの特定

時間 $$T$$ をいくつにすればコストが最小になるでしょうか。$$T=1$$ など時間が短すぎる場合、問題は実行不可能であり、したがって解はありません。

myprob = setupproblem(1);
[solm,fvalm,eflagm,outputm] = solve(myprob);

Solving problem using coneprog.
Problem is infeasible.

時間を 1.5 にすると問題が実行可能になります。

myprob = setupproblem(1.5);
[solm,fvalm,eflagm,outputm] = solve(myprob);

Solving problem using coneprog.
Optimal solution found.

この例の最後の補助関数 tomin は、時間 T の問題を設定してから、solve を呼び出して解となるコストを計算します。tomin 上で fminbnd を呼び出し、区間 $$1.5\le T\le 10$$ における最適な時間 (できる限り低いコスト) を特定します。

[Tmin,Fmin] = fminbnd(@(T)tomin(T,false),1.5,10)

Tmin = 1.9549

Fmin = 24.6102

最適時間 Tmin の軌跡を求めます。

minprob = setupproblem(Tmin);
sol = solve(minprob);

Solving problem using coneprog.
Optimal solution found.

最小となる軌跡と加速度をプロットします。

plottrajandaccel(sol,p0,pF)

最小となる解はほぼ「バンバン」解であり、加速度は 2 つの値を除くすべての値で最大かゼロのいずれかとなります。

最小とならない軌跡のプロット

さまざまな時間の軌跡をプロットします。

figure
hold on
options = optimoptions("coneprog","Display","none");
for i = 1:10
    T = 2*i;
    prob = setupproblem(T);
    sol = solve(prob,"Options",options);
    psol = sol.p;
    plot3(psol(:,1),psol(:,2),psol(:,3),'rx',"Color",[256-25*i 20*i 25*i]/256)
end
view([18 -10])
xlabel("x")
ylabel("y")
zlabel("z")
legend("T = 2","T = 4","T = 6","T = 8","T = 10","T = 12","T = 14","T = 16","T = 18","T = 20")
hold off

最も短い時間 (2) の軌跡は、このスケールではほぼ一直線となっています。時間が中程度になると、直線からは大きくずれています。最も長い時間 (20) では、再びほぼ直線となっています。

空気抵抗の導入

モデルのダイナミクスを変更して、空気抵抗を導入します。線形空気抵抗により、時間 $$t$$ 後に速度が係数 $$\exp (-t)$$ だけ変化します。運動方程式は次のようになります。

$$v(i)=v(i-1)\exp (-t)+t(a(i-1)+g)$$

$$\begin{array}{l}p(i)=p(i-1)+t\left(\frac{v(i-1)+v(i)}{2}\right)\\ =p(i-1)+t\left(\frac{1+\exp (-t)}{2}\right)v(i-1)+t^2\frac{a(i-1)+g}{2}.\end{array}$$

air = true の場合の関数 setupproblem(T,air) の問題の定式化には、速度定数を定義する行および位置定数を定義する行の両方に exp(-t) の係数があります。

vcons(2:N,:) = v(2:N,:) == v(1:(N-1),:)*exp(-t) + t*(a(1:(N-1),:) + repmat(g,N-1,1));
pcons(2:N,:) = p(2:N,:) == p(1:(N-1),:) + t*(1+exp(-t))/2*v(1:(N-1),:) + t^2/2*(a(1:(N-1),:) + repmat(g,N-1,1));

空気抵抗を含む問題に対する最適な時間を求めます。

[Tmin2,Fmin2] = fminbnd(@(T)tomin(T,true),1.5,10)

Tmin2 = 1.9394

Fmin2 = 28.7967

最適な時間は、空気抵抗がない問題の時間よりわずかに小さいだけ (1.95 ではなく 1.94) ですが、コスト Fmin は約 17% 上昇します (24.6 ではなく 28.8)。

compartable = table([Tmin;Tmin2],[Fmin;Fmin2],'VariableNames',["Time" "Cost"],'RowNames',["Original" "Air Resistance"])

compartable=2×2 table
                       Time      Cost 
                      ______    ______

    Original          1.9549     24.61
    Air Resistance    1.9394    28.797

空気抵抗を含む場合の最適解の軌跡および加速度をプロットします。

minprob = setupproblem(Tmin2,true);
sol = solve(minprob);

Solving problem using coneprog.
Optimal solution found.

plottrajandaccel(sol,p0,pF)

空気抵抗があると、ゼロ加速度の時間が後ろのタイム ステップに移動し、短縮されます。ただし、解はほぼ「バンバン」のままです。

補助関数

次のコードは、補助関数 plottrajandaccel を作成します。

function plottrajandaccel(sol,p0,pF)
figure
psol = sol.p;
plot3(psol(:,1),psol(:,2),psol(:,3),'rx')
hold on
plot3(p0(1),p0(2),p0(3),'ks')
plot3(pF(1),pF(2),pF(3),'bo')
hold off
view([18 -10])
xlabel("x")
ylabel("y")
zlabel("z")
legend("Steps","Initial Point","Final Point")
figure
asolm = sol.a;
nasolm = sqrt(sum(asolm.^2,2));
plot(nasolm,"rx")
xlabel("Time step")
ylabel("Norm(acceleration)")
end

次のコードは、補助関数 tomin を作成します。

function F = tomin(T,air)
    if nargin == 1
        air = false;
    end
    problem = setupproblem(T,air);
    opts = optimoptions("coneprog","Display","none");
    [~,F] = solve(problem,"Options",opts);
end

参考

coneprog | solve

関連するトピック

• シミュレーションまたは常微分方程式の最適化
• 常微分方程式 (ODE) の適合
• 最適化変数を使用して ODE のパラメーターを当てはめる

×

例を開く

この例の変更されたバージョンがあります。編集された方の例を開きますか?

いいえ、変更されたバージョンを上書きします。 はい