В этом примере показано, как найти положение равновесия массово-пружинной системы, висящей от двух точек привязки. Пружины имеют кусочные линейные растяжимые силы. Система состоит из массы в двух измерениях. Масса соединяется с пружинами и . Спрингс и также соединяются, чтобы разделить точки привязки. В этом случае, продолжительность нулевой силы пружины положительная длина , и пружина генерирует силу когда расширено к длине . Проблема состоит в том, чтобы найти минимальную настройку потенциальной энергии масс, куда потенциальная энергия прибывает из силы тяжести и из протяжения нелинейных пружин. Равновесие происходит в минимальной энергетической настройке.
Этот рисунок показывает пять пружин и четыре массы, приостановленные от двух точек привязки.
Потенциальная энергия массы на высоте , где ускорение свободного падения на Земле. Кроме того, потенциальная энергия идеальной линейной пружины с коэффициентом упругости расширенный к длине . Текущая модель - то, что пружина не идеальна, но имеет ненулевую продолжительность отдыха .
Математический базис этого примера прибывает от Лобо, Vandenberg, Бойда и Лебрета [1]. Для основанной на проблеме версии этого примера смотрите, Минимизируют энергию Кусочной Линейной Системы массового Spring Используя Коническое Программирование, Основанное на проблеме.
Местоположение массы , с горизонтальной координатой и вертикальная координата . Масса имеет потенциальную энергию из-за серьезности . Потенциальная энергия пружиной , где продолжительность пружины между массой и масса . Возьмите точку привязки 1 в качестве положения массы 0 и точки привязки 2 как положение массы . Предыдущее энергетическое вычисление показывает что потенциальная энергия пружины
.
При переформулировании этой проблемы потенциальной энергии, когда коническая проблема второго порядка требует введения некоторых новых переменных, как описано в Лобо [1]. Создание переменных равняйтесь квадратному корню из термина .
Пусть будьте модульным вектор-столбцом то . Проблема становится
(1)
Теперь рассмотрите как свободная векторная переменная, не данная предыдущим уравнением для . Включите отношение между и в новом наборе конических ограничений
(2)
Целевая функция еще не линейна в своих переменных, как требуется для coneprog
. Введите новую скалярную переменную . Заметьте, что неравенство эквивалентно неравенству
. (3)
Теперь проблема состоит в том, чтобы минимизировать
(4)
подвергните коническим ограничениям на и перечисленный в (2) и дополнительное коническое ограничение (3). Коническое ограничение (3) гарантирует это . Поэтому проблема (4) эквивалентна проблеме (1).
Целевая функция и конические ограничения в проблеме (4) подходят для решения с coneprog
.
Задайте шесть коэффициентов упругости , шесть констант длины , и пять масс .
k = 40*(1:6); l = [1 1/2 1 2 1 1/2]; m = [2 1 3 2 1];
Задайте аппроксимированное ускорение свободного падения на Земле .
g = 9.807;
Переменные для оптимизации являются десятью компонентами векторы, шесть компонентов вектор, и переменная. Позвольте v
будьте вектором, содержащим все эти переменные.
[v(1),v(2)]
соответствует 2D переменной .
[v(3),v(4)]
соответствует 2D переменной .
[v(5),v(6)]
соответствует 2D переменной .
[v(7),v(8)]
соответствует 2D переменной .
[v(9),v(10)]
соответствует 2D переменной .
[v(11):v(16)]
соответствует 6-D вектору .
v(17)
соответствует скалярной переменной .
Используя эти переменные, создайте соответствующий вектор целевой функции f
.
f = zeros(size(m)); f = [f;g*m]; f = f(:); f = [f;zeros(length(k)+1,1)]; f(end) = 1;
Создайте конические ограничения, соответствующие пружинам между массами (2)
.
coneprog
решатель использует конические ограничения для переменного вектора в форме
.
В следующем коде, Asc
матрица представляет термин , с bsc
= [0;0] . Коническая переменная
dsc
= и соответствующий gamma
=
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) соответствие переменная
путем создания матричного Asc
который, когда умножено на v
вектор, дает вектор . bsc
вектор соответствует постоянному 1 в термине . dsc
вектор, когда умножено на v
Возвращается . И gamma
= .
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);
Наконец, создайте нижние границы, соответствующие и переменные.
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] Лобо, Мигель Соуза, Lieven Vandenberghe, Стивен Бойд и Эрве Лебре. “Приложения Программирования Конуса Второго порядка”. Линейная алгебра и Ее Приложения 284, № 1-3 (ноябрь 1998): 193–228. https://doi.org/10.1016/S0024-3795(98)10032-0
.