このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。
錐計画法を使用した区分線形ばね質量系のエネルギーの最小化 (ソルバーベース)
この例では、2 つの固定点からぶら下がっているばね質量系の平衡位置を見つける方法を示します。ばねには区分線形張力がかかります。系は、2 つの次元の 個の質量で構成されます。質量 はばね とばね に接続されます。ばね とばね は別々の固定点に接続されます。この場合、ばね の 0 力長さは正の長さ で、ばねは長さ に伸びたときに力 を発生させます。問題は、質量の最小位置エネルギー構成を見つけることです。ここで、位置エネルギーは、重力と非線形ばねの伸びから生じます。平衡は、最小エネルギー構成で実現します。
この図は、2 つの固定点からつるされている 5 つのばねと 4 つの質量を示しています。
高さ における質量 の位置エネルギーは です。ここで、 は地球上の重力定数です。また、長さ に伸びた、ばね定数 の理想的な線形ばねの位置エネルギーは です。現在のモデルは、ばねが理想的ではありませんが、非ゼロの静止長 をもちます。
この例の数学的基礎は、Lobo、Vandenberg、Boyd、および Lebret [1] にあります。この例の問題ベースのバージョンについては、錐計画法を使用した区分線形ばね質量系のエネルギーの最小化、問題ベースを参照してください。
数学的定式化
質量 の位置は、水平座標が で、垂直座標が の です。質量 は、 の重力による位置エネルギーをもちます。ばね の位置エネルギーは です。ここで、 は、質量 と質量 間のばねの長さです。固定点 1 を質量 0 の位置として、固定点 2 を質量 の位置としてとります。前述のエネルギー計算は、ばね の位置エネルギーが以下であることを示しています。
.
この位置エネルギー問題を 2 次錐問題として再定式化するには、Lobo [1] で説明されているように、いくつかの新しい変数の導入が必要です。項 の平方根に等しい変数 を作成します。
を単位列ベクトル とします。すると、 になります。問題は次のようになります。
(1)
ここで、 を、 の以前の方程式で与えるのではなく、自由ベクトル変数と見なします。 と の関係を次の錐制約の新しいセットに組み込みます。
(2)
目的関数は、まだ、coneprog
に規定されているように、変数内で線形ではありません。新しいスカラー変数 を導入します。不等式 は次の不等式と同等であることに注意してください。
.(3)
ここで、問題は、次を最小化することです。
(4)
これは、(2) に記載されている と に対する錐制約と追加の錐制約 (3) の影響を受けます。錐制約 (3) は であることを確保します。したがって、問題 (4) は問題 (1) と同等です。
問題 (4) の目的関数と錐制約は coneprog
を使用した解法に適しています。
MATLAB® での定式化
6 つのばね制約 、6 つの長さ制約 、および 5 つの質量 を定義します。
k = 40*(1:6); l = [1 1/2 1 2 1 1/2]; m = [2 1 3 2 1];
地球 に対する近似重力定数を定義します。
g = 9.807;
最適化用の変数は、 ベクトルの 10 個の成分、 ベクトルの 6 つの成分、および 変数です。v
をこのすべての変数を含むベクトルとします。
[v(1),v(2)]
は 2 次元変数 に対応します。[v(3),v(4)]
は 2 次元変数 に対応します。[v(5),v(6)]
は 2 次元変数 に対応します。[v(7),v(8)]
は 2 次元変数 に対応します。[v(9),v(10)]
は 2 次元変数 に対応します。[v(11):v(16)]
は 6 次元ベクトル に対応します。v(17)
はスカラー変数 に対応します。
これらの変数を使用して、対応する目的関数ベクトル f
を作成します。
f = zeros(size(m)); f = [f;g*m]; f = f(:); f = [f;zeros(length(k)+1,1)]; f(end) = 1;
次の質量 (2) 間のばねに相当する錐制約を作成します。
.
coneprog
ソルバーは、次の形式の変数ベクトル の錐制約を使用します。
.
次のコードでは、Asc
行列が bsc
= [0;0]
の項 を表しています。錐変数 dsc
= と対応する gamma
=
d = zeros(1,length(f)); Asc = d; Asc([1 3]) = [1 -1]; A2 = circshift(Asc,1); Asc = [Asc;A2]; ml = length(m); dbase = 2*ml; bsc = [0;0]; for i = 2:ml gamma = -l(i); dsc = d; dsc(dbase + i) = sqrt(2/k(i)); conecons(i) = secondordercone(Asc,bsc,dsc,gamma); Asc = circshift(Asc,2,2); end
前述のコードで示したように、先端質量の位置に固定点を使用することによって、先端質量と固定点の間のばねに相当する錐制約を作成します。
x0 = [0;5]; xn = [5;4]; Asc = zeros(size(Asc)); Asc(1,(dbase-1)) = 1; Asc(2,dbase) = 1; bsc = xn; gamma = -l(ml); dsc = d; dsc(dbase + ml) = sqrt(2/k(ml)); conecons(ml + 1) = secondordercone(Asc,bsc,dsc,gamma); Asc = zeros(size(Asc)); Asc(1,1) = 1; Asc(2,2) = 1; bsc = x0; gamma = -l(1); dsc = d; dsc(dbase + 1) = sqrt(2/k(1)); conecons(1) = secondordercone(Asc,bsc,dsc,gamma);
次の 変数に相当する錐制約 (3) を作成します。
これは、行列 Asc
を作成することによって行います。この行列は、v
ベクトルを乗算すると、ベクトル が得られます。bsc
ベクトルは、項 内の定数 1 に対応します。dsc
ベクトルは、v
を乗算すると、 を返します。また、gamma
= です。
Asc = 2*eye(length(f)); Asc(1:dbase,:) = []; Asc(end,end) = -1; bsc = zeros(size(Asc,1),1); bsc(end) = -1; dsc = d; dsc(end) = 1; gamma = -1; conecons(ml+2) = secondordercone(Asc,bsc,dsc,gamma);
最後に、 変数と 変数に対応する下限を作成します。
lb = -inf(size(f)); lb(dbase+1:end) = 0;
問題の解決と解のプロット
問題の定式化は完了です。coneprog
を呼び出して問題を解きます。
[v,fval,exitflag,output] = coneprog(f,conecons,[],[],[],[],lb);
Optimal solution found.
解の点とアンカーをプロットします。
pp = v(1:2*ml); pp = reshape(pp,2,[]); pp = pp'; plot(pp(:,1),pp(:,2),'ro') hold on xx = [x0,xn]'; plot(xx(:,1),xx(:,2),'ks') xlim([x0(1)-0.5,xn(1)+0.5]) ylim([min(pp(:,2))-0.5,max(x0(2),xn(2))+0.5]) xxx = [x0';pp;xn']; plot(xxx(:,1),xxx(:,2),'b--') legend('Calculated points','Anchor points','Springs','Location',"best") hold off
パラメーターの m
、l
、および k
の値を変更して、それらがどのように解に影響するかを確認します。質量の数を変更することもできます。コードは、入力されたデータから質量の数を抽出します。
参考文献
[1] Lobo, Miguel Sousa, Lieven Vandenberghe, Stephen Boyd, and Hervé Lebret.“Applications of Second-Order Cone Programming.”Linear Algebra and Its Applications 284, no. 1–3 (November 1998):193–228. https://doi.org/10.1016/S0024-3795(98)10032-0
.