範囲制約付き二次計画法、問題ベース

ライブスクリプトを開く

この例では、二次最適化問題を解いて、サーカスのテントの形状を決定する方法を説明します。テントは、弾性のある重い素材で形成されており、制約の下で位置エネルギーが最小になる形状になります。問題を離散化すると、範囲制約付きの二次計画法問題になります。

この例のソルバーベースのバージョンについては、範囲制約付きの二次計画法、ソルバーベースを参照してください。

問題の定義

正方形の区画にサーカスのテントを建築することを考えます。このテントには 5 本の柱があり、弾力性がある重い素材で覆われます。問題は、テントの自然な形状を求めることです。位置 p におけるテントの高さ x(p) として形状をモデル化します。

高さ x まで持ち上げた重い素材の位置エネルギーは、素材の重さに比例する定数 c に対して cx です。この問題では c = 1/3000 を選択します。

素材の一部分の弾性位置エネルギー $E_{s t r e t c h}$ は、素材の高さの 2 階微分に高さを乗算した値にほぼ比例します。2 階微分は、5 つの点の有限差分近似によって近似できます (有限差分ステップのサイズは 1 であると仮定します)。 $Δ x$ は 1 番目の座標方向における 1 単位のシフトを、 $Δ y$ は 2 番目の座標方向における 1 単位のシフトを表すとします。

$E_{s t r e t c h} (p) = (- 1 (x (p + Δ_{x}) + x (p - Δ_{x}) + x (p + Δ_{y}) + x (p - Δ_{y})) + 4 x (p)) x (p) .$

テントの自然な形状では、総位置エネルギーが最小になります。問題を離散化することにより、最小化対象の総位置エネルギーは、 $E_{s t r e t c h} (p)$ + cx(p) のすべての位置 p の合計であることがわかります。

この位置エネルギーは、変数 x による二次式です。

端におけるテントの高さがゼロであるという境界条件を指定します。テントの柱の断面積は 1 x 1 単位であり、テント全体のサイズは 33 x 33 単位です。各柱の高さと位置を指定します。正方形の区画の領域とテントの柱をプロットします。

height = zeros(33);
height(6:7,6:7) = 0.3;
height(26:27,26:27) = 0.3;
height(6:7,26:27) = 0.3;
height(26:27,6:7) = 0.3;
height(16:17,16:17) = 0.5;
colormap(gray);
surfl(height)
axis tight
view([-20,30]);
title('Tent Poles and Region to Cover')

最適化問題の定式化

素材の高さを表す最適化変数 x を作成します。

x = optimvar('x',size(height));

正方形の領域の境界で x をゼロに設定します。

boundary = false(size(height));
boundary([1,33],:) = true;
boundary(:,[1,33]) = true;
x.LowerBound(boundary) = 0;
x.UpperBound(boundary) = 0;

各点の弾性位置エネルギーを計算します。はじめに、領域の内部の位置エネルギーを計算します。ここでは、解が含まれている領域を有限差分が超えることはありません。

L = size(height,1);
peStretch = optimexpr(L,L); % This initializes peStretch to zeros(L,L)
interior = 2:(L-1);
peStretch(interior,interior) = (-1*(x(interior - 1,interior) + x(interior + 1,interior) ...
    + x(interior,interior - 1) + x(interior,interior + 1)) + 4*x(interior,interior))...
    .*x(interior, interior);

領域の端では解が 0 に制約されるので、残りの項を含める必要はありません。すべての項には x の倍数が含まれており、端における x はゼロです。以下は、異なる境界条件を使用する場合の参考にするための、位置エネルギーをコメントアウトしたバージョンです。

% peStretch(1,interior) = (-1*(x(1,interior - 1) + x(1,interior + 1) + x(2,interior))...
%     + 4*x(1,interior)).*x(1,interior);
% peStretch(L,interior) = (-1*(x(L,interior - 1) + x(L,interior + 1) + x(L-1,interior))...
%     + 4*x(L,interior)).*x(L,interior);
% peStretch(interior,1) = (-1*(x(interior - 1,1) + x(interior + 1,1) + x(interior,2))...
%     + 4*x(interior,1)).*x(interior,1);
% peStretch(interior,L) = (-1*(x(interior - 1,L) + x(interior + 1,L) + x(interior,L-1))...
%     + 4*x(interior,L)).*x(interior,L);
% peStretch(1,1) = (-1*(x(2,1) + x(1,2)) + 4*x(1,1)).*x(1,1);
% peStretch(1,L) = (-1*(x(2,L) + x(1,L-1)) + 4*x(1,L)).*x(1,L);
% peStretch(L,1) = (-1*(x(L,2) + x(L-1,1)) + 4*x(L,1)).*x(L,1);
% peStretch(L,L) = (-1*(x(L-1,L) + x(L,L-1)) + 4*x(L,L)).*x(L,L);

素材の高さによる位置エネルギーを定義します。これは x/3000 です。

peHeight = x/3000;

tentproblem という名前の最適化問題を作成します。目的関数の式を含めます。これは、すべての位置における 2 つの位置エネルギーの合計です。

tentproblem = optimproblem('Objective',sum(sum(peStretch + peHeight)));

制約の設定

解が行列 height の値より大きくなければならない、という制約を設定します。この行列は、ほとんどの位置で地面を表すゼロであり、その位置におけるテントの各柱の高さが含まれています。

htcons = x >= height;
tentproblem.Constraints.htcons = htcons;

最適化ソルバーの実行

問題を解きます。出力される「ヘッシアンが対称ではありません。」という文は無視します。問題形式から二次行列への内部変換では行列が対称であることは確実ではないので、solve はこの文を出力します。

sol = solve(tentproblem);

Solving problem using quadprog.
Your Hessian is not symmetric. Resetting H=(H+H')/2.

Minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

解のプロット

最適化ソルバーで求めた解をプロットします。

surfl(sol.x);
axis tight;
view([-20,30]);