Maple を使った汎関数の最適問題

このワークシートでは、Maple を使い、汎関数(関数の関数)が最大最小となるような関数を求める問題の例をご紹介します。

例えば、

例1:
汎関数 J = 1/2; int(y(x)^2+diff(y(x),x)^2,x = 0 .. xf); が最小となるような関数 y(x) の値を求める。
ここで J を評価関数と呼び、J が最大/最小となるように y(x) を求めていきます。これを変分法を使い、解いていきます。

>    restart;

Mapleでの微分

Maple は diff コマンドや D 微分演算子を使い、微分(偏微分)を行います。例えば、次の2つは同じ偏微分を表します。

>    diff(f(x,y),x,y);

diff(f(x,y),x,y)

>    D[1,2](f)(x,y);

D[1,2](f)(x,y)

Maple では、関数で微分する (例えば、J(y(t)) を y(t) で微分する) コマンドが用意されていません。この例では多く使うので、プロシージャを用意しておくと便利です。

>    fdiff:=proc(expr,f)
algsubs(f=dummy,expr);
diff(%,dummy);
subs(dummy=f,%);
end:

>    fdiff(J(t,y(t)),y(t));

diff(J(t,y(t)),y(t))

変分法(Euler方程式の導出)

変分法は、汎関数の最大最小問題を微分方程式に置きかえる方法です。与えられた汎関数を Eulerの方程式に代入し、汎関数の最大最小問題を求めることができます。

Euler方程式の導出

積分される関数 F(x, y(x), dy(x)) とする。

>    J=Int(F(x, y(x), dy(x)),x=a..b);

J = Int(F(x,y(x),dy(x)),x = a .. b)

真の解を ty(x) と置くと、 y(x) = ty(x)+epsilon*eta(x) 。(ここで、両端条件より n(a)=n(b)=0 )。 F(x,y(x),dy(x))は?

補足: ここで求めているのは第1変分までです。近似の次数を1つ増やすと、第2変分まで求められます。ちなみに第2変分は求まる微分方程式が最大か?最小か?を決めます。

>    taylor(F(x,ty(x)+e*n(x),diff(ty(x),x)+e*diff(n(x),x)),e=0,2);

series(F(x,ty(x),D(ty)(x))+(D[2](F)(x,ty(x),D(ty)(x))*n(x)+D[3](F)(x,ty(x),D(ty)(x))*D(n)(x))*e+O(e^2),e,2)

これを使って、汎関数を書きかえると、

>    J(y(x))=J(ty(x))+epsilon*int(diff(F(x,y,dy),y)*eta(x)+diff(F(x,y,dy),dy)*diff(eta(x),x),x=a..b);

J(y(x)) = J(ty(x))+epsilon*int(diff(F(x,y,dy),y)*eta(x)+diff(F(x,y,dy),dy)*diff(eta(x),x),x = a .. b)

積分の第二項を部分積分すると、

>    with(student,intparts):

>    intparts(Int(diff(F(x,y,dy),dy)*diff(eta(x),x),x=a..b),diff(F(x,y,dy),dy));

diff(F(b,y,dy),dy)*eta(b)-diff(F(a,y,dy),dy)*eta(a)-Int(diff(F(x,y,dy),dy,x)*eta(x),x = a .. b)

ここで、第3項の x での微分は、全微分であることに注意!

両端条件より第1項と第二項は消えるので、結局、汎関数は、

>    L(y(x))=L(ty(x))+epsilon*int(collect(diff(F(x,y,dy),y)*eta(x)-diff(F(x,y,dy),dy,x)*eta(x),eta(x)),x=a..b);

L(y(x)) = L(ty(x))+epsilon*int((diff(F(x,y,dy),y)-diff(F(x,y,dy),dy,x))*eta(x),x = a .. b)

第二項がゼロであるためには、積分が常にゼロより

>    diff(F(x,y,dy),y)-`d/dx`*diff(F(x,y,dy),dy)=0;

diff(F(x,y,dy),y)-`d/dx`*diff(F(x,y,dy),dy) = 0

これが Euler の方程式!

ちなみに、F が高次の導関数を含む場合は、

F*y-diff(Fdy,x)+diff(Fddy,x,x)  - ... + (-1)^n*diff(Fdny,`$`(x,n))  = 0

Euler方程式:

例1を解いてみる

例1:

汎関数 J = 1/2   int(y(x)^2+diff(y(x),x)^2,x = 0 .. xf)  が最小となるような関数 y(x) の値を求める。

>    J=1/2*Int(y(x)^2+diff(y(x),x)^2,x = 0 .. xf);

J = 1/2*Int(y(x)^2+diff(y(x),x)^2,x = 0 .. xf)

>    F:=1/2*(y(x)^2+diff(y(x),x)^2);

F := 1/2*y(x)^2+1/2*diff(y(x),x)^2

Euler方程式を求めます。

>    deq:=fdiff(F,y(x))-diff(fdiff(F,diff(y(x),x)),x)=0;

deq := y(x)-diff(y(x),`$`(x,2)) = 0

初期値などを考慮し、求まった微分方程式から解 y(t) を求めます。

固定端点の場合: y(0)=1, y(1)=a

>    kai:=rhs(dsolve({deq,y(0)=1,y(1)=a},{y(x)}));

kai := (a-exp(-1))/(exp(1)-exp(-1))*exp(x)-1/(exp(1)-exp(-1))*(a-exp(1))*exp(-x)

解のプロット

>    a_list:=[0.9,0.4,0,-0.4]:

>    for i from 1 to nops(a_list) do
  subs(a=a_list[i],kai);
  p||i:=plot(%,x=0..1):
end do:

>    plots[display](seq(p||i,i=1..nops(a_list)),axes=box);

[Maple Plot]

自由端点:

>    kai:=rhs(dsolve({deq,y(0)=1,D(y)(1)=0},{y(x)}));

kai := exp(-1)/(exp(1)+exp(-1))*exp(x)+1/(exp(1)+exp(-1))*exp(1)*exp(-x)

>    p:=plot(kai,x=0..1,color=black):

>    plots[display]({seq(p||i,i=1..nops(a_list)),p},axes=box);

[Maple Plot]

自動的に Euler 方程式を導くプロシージャ/ deuler1, deuler2

例1のように、評価関数を Euler 方程式に当てはめていくと、この最大/最小問題を解くことができます。

これをより簡単に求められるように、プロシージャ deuler1 と deuler2 を用意しました。このプロシージャは、高階の場合と従属変数が複数の場合に対応しています。(これ以降の例は、このプロシージャを使用しています。

deuler1: deuler1( F, t, {x,y});

deuler2: deuler1( F, t, {x,y});

deuler1 と deuler2 は同じですが、deuler2 は Euler 方程式を求めるときに、1次の特別な場合を考慮します。duler1 で得られた微分方程式がうまく解けない場合に使います。deuler2で得られる微分方程式は、K1, K2,... などのパラメータが含まれるため、Mapleで直接微分方程式を解く場合には不便です。

パラメータの説明:

F:    Euler方程式に与える汎関数 F or H

t :    独立変数

{x,y}:    従属変数 (1つのときも集合で与えます。)

このプロシージャを使う前に、次のコマンドを実行し、プロシージャのファイル(deuler.m)を読み込んでください。

>    read "deuler.m";

deuler.m には、deuler1, deuler2, fdiff の3つのプロシージャが入っています。

   (注意) restart コマンドを実行すると、読み込んだプロシージャが消えてしまいます。restart コマンドを実行した後は、もう一度上のコマンドを実行してください。

いろいろな例を解いてみる

例2:

汎関数 J = int(16*y(x)^2-diff(y(x),`$`(x,2))^2+x^2,x = 0 .. xf)   が最小/最大となるような関数 y(x) の値を求める。

>    restart:
read "deuler.m";

>    F:=16*y(x)^2-diff(y(x),`$`(x,2))^2+x^2;

F := 16*y(x)^2-diff(y(x),`$`(x,2))^2+x^2

>    ode:=deuler1(F,x,{y});

ode := {32*y(x)-2*diff(y(x),`$`(x,4)) = 0}

>    dsolve(ode,y(x));

{y(x) = _C1*exp(2*x)+_C2*exp(-2*x)+_C3*sin(2*x)+_C4*cos(2*x)}

例3:

図に示される原点から直線 C=2-x にいたる曲線の長さを最小にせよ。

[Maple Plot]

>    J=Int(sqrt(1+diff(y(x),x)^2),x=0..xf);

J = Int((1+diff(y(x),x)^2)^(1/2),x = 0 .. xf)

>    F:=sqrt(1+diff(y(x),x)^2);

F := (1+diff(y(x),x)^2)^(1/2)

>    eqa:=deuler2(F,x,{y});

eqa := {y(x) = K1*x+K2}

原点を通るので、x=0 のとき y=0 より

>    solve(subs(x=0,rhs(op(eqa)))=0,{K2});

{K2 = 0}

>    eqa:=subs(%,eqa);

eqa := {y(x) = K1*x}

横断条件を考える。( F+Fy'*(C'-y')=0 )

>    C:=2-x;

C := 2-x

>    new_F:=eval(subs(eqa,F));

new_F := (1+K1^2)^(1/2)

>    F+fdiff(F,diff(y(x),x))*(diff(C,x)-diff(y(x),x))=0;

(1+diff(y(x),x)^2)^(1/2)+1/(1+diff(y(x),x)^2)^(1/2)*diff(y(x),x)*(-1-diff(y(x),x)) = 0

>    simplify(%);

-(-1+diff(y(x),x))/(1+diff(y(x),x)^2)^(1/2) = 0

>    eval(subs(eqa,-1+diff(y(x),x)=0));

-1+K1 = 0

K1=1 より、y=x

例4: 最速降下問題

粒子が点Aから出発してある曲線に沿って点Bまで動く。媒質の摩擦や抵抗を無視した場合、もっとも短い時間でAからBに到達するような曲線の形を見出せ。

[Maple Plot]

m*g*y = 1*m*v^2/2   より、

>    J=1/'sqrt(2*g)'*Int(sqrt(1+diff(y(x),x)^2)/sqrt(y(x)),x=0..xf);

J = 1/sqrt(2*g)*Int((1+diff(y(x),x)^2)^(1/2)/y(x)^(1/2),x = 0 .. xf)

>    F:=sqrt(1+diff(y(x),x)^2)/sqrt(y(x));

F := (1+diff(y(x),x)^2)^(1/2)/y(x)^(1/2)

>    ode:=deuler2(F,x,{y});

ode := {1/((1+diff(y(x),x)^2)^(1/2)*y(x)^(1/2)) = K1}

求まった微分方程式を解きます。

>    ode1:={sqrt(1+diff(y(x),x)^2)*sqrt(y(x))=K};

ode1 := {(1+diff(y(x),x)^2)^(1/2)*y(x)^(1/2) = K}

>    ode2:=y(x)=K/(1+diff(y(x),x)^2);

ode2 := y(x) = K/(1+diff(y(x),x)^2)

>    subs(diff(y(x),x)=cot(t),y=y(t),x=x(t),ode2):
ode3:=combine(convert(%,sincos));  # y(t)

ode3 := y(t)(x(t)) = 1/2*K-1/2*K*cos(2*t)

>    ode4:=diff(x(t),t)=diff(rhs(ode3),t)/cot(t);

ode4 := diff(x(t),t) = K*sin(2*t)/cot(t)

>    dsolve({ode4,x(0)=0},x(t)); # x(t)

x(t) = -1/2*K*sin(2*t)+K*t

解のプロット

>    zahyou:=[K*t-1/2*K*sin(2*t),1/2*K-1/2*K*cos(2*t)];

zahyou := [-1/2*K*sin(2*t)+K*t, 1/2*K-1/2*K*cos(2*t)]

>    K:=1;

K := 1

>    plot([zahyou[1],-zahyou[2],t=0..Pi],scaling=constrained);

[Maple Plot]

求まった曲線は、サイクロイド曲線です。

条件付き変分問題

汎関数の最大/最小問題には、求める関数にいろいろな制約条件がつくことがあります。

そのような場合は、F に制約条件を考慮した H をeuler方程式に当てはめます。

積分型制約条件( Int(G(x,y,dy),x) = l   )の場合: H = F+lambda*G  

関数型制約条件( g(x,y) = 0  )の場合: H = F+lambda(t)*g  

微分方程式型制約条件( g(x,y,dy) = 0  )の場合: H = F+lambda(t)*g  

lambda, lambda(t)  は、ラグランジュ乗数と呼ばれます。

例5: 最大面積問題

長さ L のひもを y(0)=0, y(1)=0 で固定したとき、x 軸とひもで囲む面積が最大になるような ひもの関数 y(x) を求めよ。

[Maple Plot]

>    # 評価関数
J=Int(y(x),x=0..1);

J = Int(y(x),x = 0 .. 1)

>    # 制約条件
L=Int(sqrt(1+diff(y(x),x)^2),x=0..1);

L = Int((1+diff(y(x),x)^2)^(1/2),x = 0 .. 1)

>    H:=y(x)+lambda*sqrt(1+diff(y(x),x)^2);

H := y(x)+lambda*(1+diff(y(x),x)^2)^(1/2)

>    ode:=deuler1(H,x,{y});

ode := {-(-(1+diff(y(x),x)^2)^(1/2)-(1+diff(y(x),x)^2)^(1/2)*diff(y(x),x)^2+lambda*diff(y(x),`$`(x,2)))/(1+diff(y(x),x)^2)^(3/2) = 0}

>    kai:=dsolve(ode union {y(0)=0,y(1)=0},y(x));

kai := y(x) = -(-x^2+x-1/4+lambda^2)^(1/2)+1/2*(-1+4*lambda^2)^(1/2), y(x) = (-x^2+x-1/4+lambda^2)^(1/2)-1/2*(-1+4*lambda^2)^(1/2)

>    yx:=rhs(kai[1]);

yx := -(-x^2+x-1/4+lambda^2)^(1/2)+1/2*(-1+4*lambda^2)^(1/2)

>    lambda:=0.5;

lambda := .5

>    plot([x,yx,x=0..1],scaling=constrained);

[Maple Plot]

>    lambda:='lambda':

>   

最適制御問題

制御では、あるものの振る舞いが微分方程式で与えられ、ある条件を最大/最小にする入力 u(t) を求めるという問題があります。

この問題を、最適制御問題と呼びます。

この場合、微分方程式が制約条件となり、条件付き変分問題として扱う事ができます。(上の章:微分方程式型制約条件( g(x,y,dy) = 0  )の場合: H = F+lambda(t)*g  )

例6:

評価関数:

J = int(x(t)^2+u(t)^2,t = 0 .. tf)

状態方程式:

diff(x(t),t) = -x(t)+u(t) , x(0) = 1 , tf = 1/2 , x(tf) = 0       

>    H:= (x(t)^2+u(t)^2) - lambda(t)*( diff(x(t),t)+x(t)-u(t));

H := x(t)^2+u(t)^2-lambda(t)*(diff(x(t),t)+x(t)-u(t))

>    kai:=deuler1(H,t,{x,u,lambda});

kai := {lambda(t)+2*u(t) = 0, -lambda(t)+2*x(t)+diff(lambda(t),t) = 0, -diff(x(t),t)-x(t)+u(t) = 0}

>    deqs:=map2(algsubs,u(t)=-lambda(t)/2,{-diff(x(t),t)-x(t)+u(t) = 0, -lambda(t)+2*x(t)+diff(lambda(t),t) = 0});

deqs := {-diff(x(t),t)-x(t)-1/2*lambda(t) = 0, -lambda(t)+2*x(t)+diff(lambda(t),t) = 0}

>    evalf(dsolve(deqs union {x(0)=1,x(1/2)=0} ,{x(t),lambda(t)}));

{x(t) = -.3212077023*exp(1.414213562*t)+1.321207702*exp(-1.414213562*t), lambda(t) = 1.550927983*exp(1.414213562*t)+1.094524298*exp(-1.414213562*t)}
{x(t) = -.3212077023*exp(1.414213562*t)+1.321207702*exp(-1.414213562*t), lambda(t) = 1.550927983*exp(1.414213562*t)+1.094524298*exp(-1.414213562*t)}

>    plot(subs(%,x(t)),t=0..1,axes=boxed);

[Maple Plot]

>   

ハミルトン関数

ハミルトンの関数を使い、例6の問題を解いてみます。(より一般的な解法です。求める微分方程式が1次の形で求まります。)

例7:

評価関数:

J = int(x(t)^2+u(t)^2,t = 0 .. tf)

状態方程式:

diff(x(t),t) = -x(t)+u(t) , x(0) = 1 , tf = 1/2 , x(tf) = 0       

>    H0:=x(t)^2+u(t)^2+lambda(t)*(-x(t)+u(t));

H0 := x(t)^2+u(t)^2+lambda(t)*(-x(t)+u(t))

>    solve(fdiff(H0,u(t))=0,{u(t)});

{u(t) = -1/2*lambda(t)}

>    maruH:=subs(%,H0);

maruH := x(t)^2+1/4*lambda(t)^2+lambda(t)*(-x(t)-1/2*lambda(t))

>    de1:= diff(x(t),t)=fdiff(maruH,lambda(t));

de1 := diff(x(t),t) = -x(t)-1/2*lambda(t)

>    de2 := diff(lambda(t),t)=-fdiff(maruH,x(t));

de2 := diff(lambda(t),t) = lambda(t)-2*x(t)

>    evalf(dsolve({de1,de2,x(0)=1,x(1/2)=0},{x(t),lambda(t)}));

{x(t) = -.3212077023*exp(1.414213562*t)+1.321207702*exp(-1.414213562*t), lambda(t) = 1.550927983*exp(1.414213562*t)+1.094524298*exp(-1.414213562*t)}
{x(t) = -.3212077023*exp(1.414213562*t)+1.321207702*exp(-1.414213562*t), lambda(t) = 1.550927983*exp(1.414213562*t)+1.094524298*exp(-1.414213562*t)}

>    plot(subs(%,x(t)),t=0..1,axes=boxed);

[Maple Plot]

>   

ロケットカーの最短時間問題

止まっているロケットカーを、ある距離 d 動かすときの、最短時間を求めます。

>    restart;

ロケットカーのシステムは、次のような式で表すことができます。 x(t) は車の位置、u(t) はアクセルによる入力です。 また、-1< u(t) < 1 とし、ブレーキの強さとアクセルの強さは同じであると仮定します。

>    diff(x(t),t,t)=u(t);

diff(x(t),`$`(t,2)) = u(t)

最短時間問題の評価関数です。

>    J=Int(1,t=0..T);

J = Int(1,t = 0 .. T)

制御入力 u(t) に制限があるこのような問題は、バング・バング制御と呼ばれ、入力関数は次の式を使います。(実際には、ポントリャギンの最小原理を使うことにより求まります。)

>    assume(T>0);

>    u:=piecewise(t>=0 and t<T/2, 1, t>=T/2 and t<T, -1);

u := PIECEWISE([1, 0 <= t and t < 1/2*T],[-1, 1/2*T <= t and t < T])

この入力を使い、x(t) を求めると。

>    ans:=dsolve({ diff(x(t),t,t)=u,x(0)=0,D(x)(0)=0},x(t));

ans := x(t) = PIECEWISE([0, t < 0],[1/2*t^2, t <= 1/2*T],[-1/2*t^2-1/4*T^2+T*t, t <= T],[1/4*T^2, T < t])

結果を見ると、t>Tでは定数となっていて、時間Tで車が止まります。T^2/4; は進んだ距離です。

>    solve(d=T^2/4,{T});

{T = 2*d^(1/2)}, {T = -2*d^(1/2)}

距離 d 進むための最短時間は 2*sqrt(d); です。

プロット:

>    x:=subs(ans,x(t));

x := PIECEWISE([0, t < 0],[1/2*t^2, t <= 1/2*T],[-1/2*t^2-1/4*T^2+T*t, t <= T],[1/4*T^2, T < t])

>    x:=map2(subs,T=4,x);

x := PIECEWISE([0, t < 0],[1/2*t^2, t <= 2],[-1/2*t^2-4+4*t, t <= 4],[4, 4 < t])

>    plot(x,t=0..4);

[Maple Plot]

位相のプロット:

>    diff(x,t);

PIECEWISE([0, t <= 0],[t, t <= 2],[4-t, t <= 4],[0, 4 < t])

>    plot([x,diff(x,t),t=0..4]);

[Maple Plot]

>   

ロケットカーの最小エネルギー問題

止まっているロケットカーを、ある距離 d 動かすとき、消費するエネルギーが最小となるような 制御入力 u(t) を求めます。

>    restart;

ロケットカーのシステムは、次のような式で表すことができます。 x(t) は車の位置、u(t) はアクセルとブレーキによる制御入力です。 また、-1< u(t) < 1 とし、ブレーキの強さとアクセルの強さは同じであると仮定します

>    diff(x(t),t,t)=u(t);

diff(x(t),`$`(t,2)) = u(t)

エネルギーを最小にする評価関数です。

>    J=Int(k+abs(u),t=0..T);

J = Int(k+abs(u),t = 0 .. T)

ピースワイズ関数を使い、入力 u(t) を定義します。入力 u(t) は、最大のアクセルを時間 b まで、その後 T-b まで惰性で動かし、T まで最大のブレーキをかけるという形にします。(このように u(t)
を数式の形で与えることにより、最適制御の問題は、パラメータの最適問題に変換できます。)

>    assume(b>0,b<T/2);

>    u:=piecewise(t>=0 and t<b, 1, t>T-b and t<T, -1);

u := PIECEWISE([1, 0 <= t and t < b],[-1, T-b < t and t < T])

この入力を使い、x(t) を求めると。

>    ans:=dsolve({ diff(x(t),t,t)=u,x(0)=0,D(x)(0)=0},x(t));

ans := x(t) = PIECEWISE([0, t < 0],[1/2*t^2, t < b],[-1/2*b^2+b*t, t < T-b],[-b^2-1/2*t^2-1/2*T^2+T*b+T*t, t < T],[-b^2+T*b, T <= t])

結果から 進んだ距離 d は、d = T*b-b^2  

>    'T'=solve(d=b*T-b^2,T);

T = (d+b^2)/b

コスト関数を計算すると。

>    J:=int(k+abs(u),t=0..T);

J := 2*k*b+2*b+k*(T-2*b)

d=3, k=1/8 とするとコスト関数は、

>    Tans:=solve(d=b*T-b^2,T);

Tans := (d+b^2)/b

>    subs(T=%,J): JJ:=subs(d=3,k=1/8,%);

JJ := 2*b+1/8*(3+b^2)/b

この JJ が最小になるように b を求めます。

>    solve(diff(JJ,b)=0); evalf(%);

1/17*51^(1/2), -1/17*51^(1/2)

.4200840252, -.4200840252

プロット:

>    x:=subs(ans,T=Tans,d=3,b=%[1],x(t));

x := PIECEWISE([0, t < 0],[1/2*t^2, t < .4200840252],[-.8823529410e-1+.4200840252*t, t < 7.141428428],[-1/2*t^2-25.58823530+7.561512453*t, t < 7.561512453],[3, 7.561512453 <= t])

>    plot([x,diff(x,t),t=0..8]);

[Maple Plot]

また k をパラメータのままにして、fuel efficiency factor の与える影響を見てみましょう。

>    subs(T=Tans,J): JJ:=subs(d=3,%);

JJ := 2*k*b+2*b+k*((3+b^2)/b-2*b)

この JJ が最小になるように b を求めます。

>    solve(diff(JJ,b)=0,{b});

{b = 1/(k+2)*3^(1/2)*((k+2)*k)^(1/2)}, {b = -1/(k+2)*3^(1/2)*((k+2)*k)^(1/2)}

>    x:=subs(ans,T=Tans,d=3,%[1],x(t));

x := PIECEWISE([0, t < 0],[1/2*t^2, t < 1/(k+2)*3^(1/2)*((k+2)*k)^(1/2)],[-3/2*k/(k+2)+1/(k+2)*3^(1/2)*((k+2)*k)^(1/2)*t, t < 1/3*(3+3*k/(k+2))*(k+2)*3^(1/2)/((k+2)*k)^(1/2)-1/(k+2)*3^(1/2)*((k+2)*k)^(...

>    for i from 1 to 10 do
  xx:=evalf(subs(k=i/10, x));
  p||i:=plot([xx,diff(xx,t),t=0..10]):
end do:

>    plots[display]([seq(p||i,i=1..10)],insequence=true);

[Maple Plot]

fuel efficiency factor (k)が大きくなると、アクセルとブレーキを踏む時間が多くなる。

>   

>