Математическая формулировка

Местоположение массы $$i$$ - $$x(i$$), с горизонтальной координатой $${}_{\mathrm{}}x1($$i) и вертикальной $${\mathrm{координатой}}_{\mathrm{}}\mathrm{x2}$$(i). $$$$Масса i имеет потенциальную энергию, обусловленную гравитацией $$gm{}_{\mathrm{(}}i)$$x2 (i). Потенциальная энергия $$$$в пружине $$iравнаk(i)(d{(}^{\mathrm{i}})\mathrm{}$$-1 (i)) $$\mathrm{2/2},$$ где d (i) - длина пружины $$$$между массой $$i\mathrm{}$$и массой i-1. Возьмем точку привязки 1 в качестве положения массы 0, а точку привязки 2 в качестве $$\mathrm{положения}$$массы n + 1. Предыдущий расчет энергии показывает, что потенциальная $$$$энергия пружины i равна

$$\mathrm{Энергия}(i)\frac{=k(i)(\Vert x(\mathrm{i})-x({}^{\mathrm{i-1}}}{\mathrm{)}}$$‖ -l (i)) 22.

Переформулирование этой потенциальной энергетической проблемы как проблемы программирования конуса второго порядка требует введения некоторых новых переменных, как описано в Лобо [1]. Создайте переменные $$t(i$$), равные квадратному корню термина $$Энергия($$i).

Пусть $$e$$ - вектор единичного столбца $$[\genfrac{}{}{0ex}{}{\mathrm{}}{01}$$]. Затем $${}_{\mathrm{x2}}(i){}^{=}eTx$$(i). Проблема становится

$$\underset{\mathrm{minx},}{\mathrm{}}t\underset{(}{}\mathrm{\sum igm}{(}^i)eTx(i^{\mathrm{)}}$$+‖t‖2      ).         (1)

Теперь рассмотрим $$t$$ как переменную свободного вектора, не заданную предыдущим уравнением для $$t(i$$). Включить взаимосвязь между $$x($$i) $$иt$$(i) в новый набор ограничений конуса

$$\Vert x(i)-x\mathrm{}(\mathrm{i-1})\Vert \sqrt{\frac{\mathrm{-}}{1(i}})\le 2k$$  ( i) t (i). (2)

Целевая функция еще не является линейной в своих переменных, как требуется для coneprog. Введите новую скалярную переменную $$y$$. Обратите внимание, что неравенство $$\Vert {}^{\mathrm{}}$$t‖2≤y эквивалентно неравенству

$$\Vert \genfrac{}{}{0ex}{}{\mathrm{[}}{\mathrm{}}2t1-y]\mathrm{\Vert }$$≤1+y   .                            (3)

Теперь проблема сводится к минимуму

$$\underset{\mathrm{minx},t,}{\mathrm{}}\underset{}{y}(\sum igm{}^{(}i)\mathrm{eTx}$$ ( i )  +  y )  ( 4 )         

с учетом ограничений конуса на $$x(i$$) и $$t($$i), перечисленных в пункте (2), и дополнительного ограничения конуса (3). Зависимость конуса (3) обеспечивает $${\Vert }^{\mathrm{}}$$t‖2≤y. Следовательно, проблема (4) эквивалентна проблеме (1).

Целевая функция и ограничения конуса в задаче (4) подходят для решения с coneprog.

Состав MATLAB ®

Определите шесть пружинных констант $$k$$, шесть констант длины $$l$$и пять масс $$м$$.

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

Определите переменные оптимизации, соответствующие математическим переменным задачи. Для простоты задайте опорные точки как две виртуальные массовые точки x(1,:) и x(end,:). Эта композиция позволяет каждой пружине растягиваться между двумя массами.

nmass = length(m) + 2;
% k and l have nmass-1 elements
% m has nmass - 2 elements
x = optimvar('x',[nmass,2]);
t = optimvar('t',nmass-1,'LowerBound',0);
y = optimvar('y','LowerBound',0);

Создайте задачу оптимизации и задайте для целевой функции выражение в (4).

prob = optimproblem;
obj = dot(x(2:(end-1),2),m)*g + y;
prob.Objective = obj;

Создайте зависимости конуса, соответствующие выражению (2).

conecons = optimineq(nmass - 1);
for ii = 1:(nmass-1)
    conecons(ii) = norm(x(ii+1,:) - x(ii,:)) - l(ii) <= sqrt(2/k(ii))*t(ii);
end
prob.Constraints.conecons = conecons;

Укажите опорные точки anchor0 и anchorn. Создайте ограничения равенства, указывающие, что две виртуальные конечные массы расположены в опорных точках.

anchor0 = [0 5];
anchorn = [5 4];
anchorcons = optimeq(2,2);
anchorcons(1,:) = x(1,:) == anchor0;
anchorcons(2,:) = x(end,:) == anchorn;
prob.Constraints.anchorcons = anchorcons;

Создайте зависимость конуса, соответствующую выражению (3).

ycone = norm([2*t;(1-y)]) <= 1 + y;
prob.Constraints.ycone = ycone;

Решить проблему

Постановка проблемы завершена. Решить проблему, позвонив solve.

[sol,fval,eflag,output] = solve(prob);

Solving problem using coneprog.
Optimal solution found.

Постройте график точек решения и анкеров.

plot(sol.x(2:(nmass-1),1),sol.x(2:(nmass-1),2),'ro')
hold on
plot([sol.x(1,1),sol.x(end,1)],[sol.x(1,2),sol.x(end,2)],'ks')
plot(sol.x(:,1),sol.x(:,2),'b--')
legend('Calculated points','Anchor points','Springs','Location',"best")
xlim([sol.x(1,1)-0.5,sol.x(end,1)+0.5])
ylim([min(sol.x(:,2))-0.5,max(sol.x(:,2))+0.5])
hold off

Можно изменить значения параметров m, l, и k чтобы увидеть, как они влияют на решение. Можно также изменить количество формообразующих элементов; код берет количество масс из предоставленных данных.

Ссылки

[1] Лобо, Мигель Соуза, Ливен Ванденберге, Стивен Бойд и Эрве Лебрет. «Приложения программирования конуса второго порядка». Линейная алгебра и ее приложения 284, № 1-3 (ноябрь 1998): 193-228. https://doi.org/10.1016/S0024-3795(98)10032-0.