这个例子展示了如何通过求解一个二次优化问题来确定马戏团帐篷的形状。帐篷由沉重的弹性材料制成,并形成一个受约束的势能最小的形状。问题的离散化可得到一个有界约束的二次规划问题。
有关此示例的基于求解器的版本,请参见基于求解器的有界约束二次规划.
考虑建造马戏团帐篷覆盖平方。帐篷有五个焊接覆盖着沉重的弹性材料。问题是找到帐篷的自然形状。将形状模拟为高度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初始化为零(L,L)室内= 2:(l - 1);peStretch = (-1*(x(interior -1,interior) + x(interior + 1,interior))...+ x(interior,interior - 1) + x(interior,interior + 1)) + 4*x(interior,interior)).... * 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,室内);= (-1*(x(interior - 1,1) + x(interior + 1,1) + x(interior,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 =解决(帐篷问题);
用四轴飞行器解决问题。你的Hessian不是对称的。重置H = (H + H) / 2。找到满足约束条件的最小值。优化完成是因为目标函数在可行方向上不减小到最优性公差的值内,约束条件满足到约束公差的值内。
绘制优化求解器发现的解决方案。
surfl (sol.x);轴紧的;查看([ - 20,30]);