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

Местоположение массы $$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$$.

g = 9.807;

Переменные для оптимизации являются десятью компонентами $$x$$ векторы, шесть компонентов $$t$$ вектор, и $$y$$ переменная. Позвольте v будьте вектором, содержащим все эти переменные.

Используя эти переменные, создайте соответствующий вектор целевой функции f.

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

Создайте конические ограничения, соответствующие пружинам между массами (2)

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

coneprog решатель использует конические ограничения для переменного вектора $$v$$ в форме

$$\Vert A_{sc}\cdot v-b_{sc}\Vert \le d_{sc}^{T}v-\gamma$$.

В следующем коде, Asc матрица представляет термин $$\Vert x(i)-x(i-1)\Vert$$, с bsc= [0;0] . Коническая переменная dsc = $$\sqrt{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);

Создайте коническое ограничение (3) соответствие $$y$$ переменная

путем создания матричного Asc который, когда умножено на v вектор, дает вектор $$\left[\genfrac{}{}{0ex}{}{2t}{-y}\right]$$. bsc вектор соответствует постоянному 1 в термине $$1-y$$. 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);

Наконец, создайте нижние границы, соответствующие $$t$$ и $$y$$ переменные.

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

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

Ссылки

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