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

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

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\mathrm{}(\mathrm{i-1})\Vert \sqrt{\frac{\mathrm{-}}{1(i}})$$≤2k (i) t (i).

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

$$\Vert {}_{}{}_{}{}_{}^{}$$Asc⋅v-bsc‖≤dscTv-γ.

В следующем коде Asc матрица представляет член $$\Vert x(i)-x\mathrm{}($$i-1) ‖, сbsc = [0;0]. Конусная переменная dsc = $$\sqrt{\mathrm{}\mathrm{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 вектор, дает вектор $$[\genfrac{}{}{0ex}{}{\mathrm{}}{\mathrm{2t-y}}$$]. bsc вектор соответствует константе 1 в члене $$\mathrm{}\mathrm{1-y}$$. dsc вектор, при умножении на v, возвращает $$y$$. И gamma $$\mathrm{}$$= -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

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

Ссылки

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