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

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

$$Energy(i)=\frac{k(i)(\Vert x(i)-x(i-1)\Vert -l(i){)}^2}{2}$$.

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

Пусть $$e$$ будьте модульным вектор-столбцом $$\left[\genfrac{}{}{0ex}{}{0}{1}\right]$$то $$x_2(i)=e^{T}x(i)$$. Проблема становится

$$\underset{x,t}{\min }(\sum _{i}gm(i)e^{T}x(i)+\Vert t{\Vert }^2).$$               (1)

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

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

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

$$\Vert \left[\genfrac{}{}{0ex}{}{2t}{1-y}\right]\Vert \le 1+y$$                               . (3)

Теперь проблема состоит в том, чтобы минимизировать

$$\underset{x,t,y}{\min }(\sum _{i}gm(i)e^{T}x(i)+y)$$                     (4)

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

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

MATLAB® Formulation

Задайте шесть коэффициентов упругости $$k$$, шесть констант длины $$l$$, и пять масс $$m$$.

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

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

Ссылки

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