这个例子展示了如何通过解决二次优化问题来确定马戏团帐篷的形状。帐篷由重的、有弹性的材料制成,并在限制条件下形成势能最小的形状。将问题离散化,得到有界约束二次规划问题。
有关此示例的基于求解器的版本,请参见基于求解器的有界约束二次规划.
考虑建造马戏团帐篷覆盖平方。帐篷有五个焊接覆盖着沉重的弹性材料。问题是找到帐篷的自然形状。将形状模拟为高度x(p)帐篷在位置p.
重物上升到高处的势能x是残雪,常数c这与材料的重量成比例。对于这个问题,选择c= 1/3000。
材料的弹性势能 近似正比于物质高度的二阶导数,乘以高度。你可以用五点有限差分近似来近似二阶导数(假设有限差分步骤的大小为1) 表示在第一坐标方向上1的偏移, 表示在第二个坐标方向上的移动1。
帐篷的自然形状最小化总势能。通过离散解决问题,您发现最小化的总潜在能量是所有位置的总和p的 +残雪(p)。
这种势能是变量中的二次表达x
.
指定边界条件,即帐篷在边缘处的高度为零。帐篷柱子的横截面是1 × 1的单元,帐篷的总尺寸是33 × 33的单元。指定每根杆的高度和位置。绘制方形地块区域和帐篷柱。
身高= 0 (33);高度(者者)= 0.3;高度(二六27,二六27)= 0.3;高度(者,二六27)= 0.3;高度(二六27,者)= 0.3;高度(16:17,16:17)= 0.5;colormap(灰色);surfl(高度)轴紧的查看([ - 20,30]);标题(“帐篷支柱和需要覆盖的区域”)
创建优化变量x
表示材料的高度。
x = optimvar ('X',尺寸(高度));
集x
在方框的边界上为零。
边界= false(大小(高度));边界([1,33]:)= true;边界(:,(1,33))= true;x.LowerBound(边界)= 0;x.UpperBound(边界)= 0;
计算每一点的弹性势能。首先,计算区域内部的势能,其中有限差分不超过包含解的区域。
L =大小(高度,1);peStretch = optimexpr (L, L);初始化peStretch为0 (L,L)室内= 2:(l - 1);peStretch(inner, inner) = (-1*(x(inner -1, inner) + x(inner + 1, inner))...+ x(内部,内部- 1)+ x(内部,内部+ 1)+ 4*x(内部,内部).... * x(室内,室内);
因为解在区域的边缘被约束为0,你不需要包含剩下的项。所有项都有倍数x
, 和x
在边缘为零。有关您想要使用不同的边界条件的情况下,以下是潜在能量的注释版本。
%Pertrach(1,内部)=(-1 *(x(1,内部 - 1)+ x(1,内部+ 1)+ x(2,内部))...% + 4 * x(内部))。* x(内部);%Pertrach(l,内部)=(-1 *(x(l,内部 - 1)+ x(l,内部+ 1)+ x(l-1,内部))...% + 4 * x (L,内政部))* x (L,室内);% peStretch(inner,1) = (-1*(x(inner - 1,1) + x(inner + 1,1) + x(inner,2))…% + 4 * x(内部,1))* x(内部,1);%Pertrch(内部,L)=(-1 *(x(内部 - 1,l)+ x(内部+ 1,l)+ x(内部,l-1))...% + 4 * x(内部,L)) * x(室内,L);% peStretch (1, - 1) = (1 * (x (2, 1) + x(1、2))+ 4 * x(1,1))。* x (1,1);% peStretch (L) = (1 * (x (2 L) + x (1, L - 1)) + 4 * x(1升))。* x(1升);% 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;
创建一个名为帐篷问题
.包括目标函数的表达式,这是所有位置的两个潜在能量的总和。
TentProblus = OptimProblem(“目标”总和(sum (peStretch + peHeight)));
设置解必须位于值之上的约束高度
矩阵。该矩阵在大多数位置处为零,代表地面,并且包括其位置处的每个帐篷杆的高度。
Htcons = x >= height;tentproblem.Constraints.htcons = htcons;
解决这个问题。忽略结果语句“Your Hessian is not symmetric.”。解决
发出此语句,因为从问题表单到二次矩阵的内部转换不会确保矩阵是对称的。
sol =解决(帐篷问题);
使用quadprog解决问题。你的黑森不是对称的。重置H = (H + H) / 2。求满足约束条件的最小值。优化完成是因为目标函数在可行方向上是不递减的,到最优公差的值以内,并且约束满足到约束公差的值以内。
绘制优化求解器发现的解决方案。
surfl (sol.x);轴紧的;查看([ - 20,30]);