基于问题的有界约束二次规划

此示例说明如何通过求解二次优化问题来确定马戏团帐篷的形状。帐篷由厚重的弹性材料制成，在约束下形成具有最小势能的形状。问题的离散化导致了有界约束的二次规划问题。

问题定义

考虑搭建一个马戏团帐篷来覆盖一块方形场地。帐篷有五根杆子，上面覆盖着厚重的弹性材料。问题在于找到帐篷的自然形状。将形状建模为位置 p 处帐篷的高度 x(p)。

重物被提升到高度x的势能为cx，其中c为常数，与物的重量成比例。对于这个问题，选择 c = 1/3000。

一块材料 $E_{s t r e t c h}$ 的弹性势能大约与材料高度的二阶导数乘以高度成正比。您可以通过五点有限差分近似来近似二阶导数（假设有限差分步骤的大小为 1）。让 $Δ x$ 表示在第一个坐标方向上平移 1，而 $Δ y$ 表示在第二个坐标方向上平移 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×1 单位，帐篷的总尺寸为 33×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 的优化问题。包括目标函数的表达式，它是所有位置的两个势能的总和。

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

设置约束

设置约束，即解必须位于 height 矩阵的值之上。该矩阵在大多数位置为零，代表地面，并包括每个帐篷杆在其位置的高度。

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

运行优化求解器

求解。忽略结果语句“您的 Hessian 不对称。” 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]);