Main Content

使用基于求解器的锥规划最小化分段线性质量-弹簧系统的能量

此示例显示如何找到悬挂在两个锚点的质量-弹簧系统的平衡位置。弹簧具有分段线性拉力。该系统由二维空间中的 n 个质量组成。质量 i 与弹簧 ii+1 相连。弹簧 1n+1 也连接到单独的锚点。在这种情况下,弹簧 i 的零力长度为正长度 l(i),并且弹簧在拉伸到长度 q+l(i) 时会产生力 k(i)q。问题在于找到质量的最小势能配置,其中势能来自重力和非线性弹簧的拉伸。平衡发生在最小能量配置上。

该图显示了五个弹簧和四个质量悬挂在两个锚点上。

质量 m 在高度 h 处的势能为 mgh,其中 g 是地球引力常数。另外,弹簧常数为 k 的理想线性弹簧拉伸至长度 q 时,其势能为 kq2/2。当前的模型是弹簧不是理想的,但具有非零的静止长度 l

本例的数学基础来自 Lobo、Vandenberg、Boyd 和 Lebret [1]。有关此示例的基于问题的版本,请参阅 基于问题的锥规划最小化分段线性质量-弹簧系统的能量

数学表示

质量 i 的位置为 x(i),水平坐标为 x1(i),垂直坐标为 x2(i)。质量 i 由于 gm(i)x2(i) 的重力而具有势能。弹簧 i 中的势能为 k(i)(d(i)-l(i))2/2,其中 d(i) 是质量 i 和质量 i-1 之间的弹簧长度。以锚点 1 作为质量 0 的位置,以锚点 2 作为质量 n+1 的位置。前面的能量计算表明,弹簧 i 的势能为

Energy(i)=k(i)(x(i)-x(i-1)-l(i))22.

将此势能问题重新表述为二阶锥问题需要引入一些新变量,如 Lobo [1] 中所述。创建等于项 Energy(i) 的平方根的变量 t(i)

t(i)=k(i)(x(i)-x(i-1)-l(i))22.

e 为单位列向量 [01]。则 x2(i)=eTx(i)。问题变成了

minx,t(igm(i)eTx(i)+t2). (1)

现在将 t 视为一个自由向量变量,不是由前面的 t(i) 方程给出的。将 x(i)t(i) 之间的关系纳入新的锥约束

x(i)-x(i-1)-l(i)2k(i)t(i). (2)

目标函数在其变量中还不是线性的,正如 coneprog 所要求的那样。引入一个新的标量变量 y。请注意,不等式 t2y 等价于不等式

[2t1-y]1+y.(3)

现在的问题是最小化

minx,t,y(igm(i)eTx(i)+y) (4)

受制于 (2) 中列出的 x(i)t(i) 的锥约束以及 (3) 的附加锥约束。锥约束 (3) 确保 t2y。因此,问题(4)等价于问题(1)。

问题 (4) 中的目标函数和锥约束适合用 coneprog 解。

MATLAB® 表示

定义六个弹簧常数 k、六个长度常数 l 和五个质量 m

k = 40*(1:6);
l = [1 1/2 1 2 1 1/2];
m = [2 1 3 2 1];

定义地球上的近似引力常数 g

g = 9.807;

需要优化的变量是 x 向量的十个分量、t 向量的六个分量以及 y 变量。让 v 成为包含所有这些变量的向量。

  • [v(1),v(2)] 对应于二维变量 x(1)

  • [v(3),v(4)] 对应于二维变量 x(2)

  • [v(5),v(6)] 对应于二维变量 x(3)

  • [v(7),v(8)] 对应于二维变量 x(4)

  • [v(9),v(10)] 对应于二维变量 x(5)

  • [v(11):v(16)] 对应于 6 维向量 t

  • v(17) 对应于标量变量 y

使用这些变量,创建相应的目标函数向量 f

f = zeros(size(m));
f = [f;g*m];
f = f(:);
f = [f;zeros(length(k)+1,1)];
f(end) = 1;

创建与质量之间的弹簧相对应的锥约束(2)

x(i)-x(i-1)-l(i)2k(i)t(i).

coneprog 求解器使用锥约束来求解变量向量 v,形式如下

Ascv-bscdscTv-γ.

在下面的代码中,Asc 矩阵代表项 x(i)-x(i-1),其中 bsc = [0;0]。锥变量 dsc = 2/k(i) 以及对应的 gamma = -l(i).

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);

创建与 y 变量对应的锥约束(3)

[2t1-y]1+y

通过创建矩阵 Asc ,将其与 v 向量相乘,得到向量 [2t-y]bsc 向量对应于项 1-y 中的常数 1。dsc 向量与 v 相乘后返回 y。并且 gamma = -1

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);

最后,创建与 ty 变量对应的下界。

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

您可以更改参数 mlk 的值,看看它们如何影响解。您还可以改变质量数;代码会从您提供的数据中获取质量数。

参考资料

[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.

另请参阅

|

相关主题