Минимизируйте энергию кусочной линейной системы массового Spring Используя коническое программирование

В этом примере показано, как найти положение равновесия массово-пружинной системы, висящей от двух точек привязки. Пружины имеют кусочные линейные растяжимые силы. Система состоит из n массы в двух измерениях. Масса i соединяется с пружинами i и i+1. Спрингс 1 и n+1 также соединяются, чтобы разделить точки привязки. В этом случае, продолжительность нулевой силы пружины i положительная длина l(i), и пружина генерирует силу k(i)q когда расширено к длине q+l(i). Проблема состоит в том, чтобы найти минимальную настройку потенциальной энергии масс, куда потенциальная энергия прибывает из силы тяжести и из протяжения нелинейных пружин. Равновесие происходит в минимальной энергетической настройке.

Этот рисунок показывает пять пружин и четыре массы, приостановленные от двух точек привязки.

Потенциальная энергия массы m на высоте h mgh, где g ускорение свободного падения на Земле. Кроме того, потенциальная энергия идеальной линейной пружины с коэффициентом упругости k расширенный к длине q kq2/2. Текущая модель - то, что пружина не идеальна, но имеет ненулевую продолжительность отдыха l.

Математическое основание этого примера прибывает от Лобо, Vandenberg, Бойда и Лебрета [1].

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

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

Energy(i)=k(i)(x(i)-x(i-1)-l(i))22.

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

t(i)=k(i)(x(i)-x(i-1)-l(i))22.

Пусть e будьте модульным вектор-столбцом [01]то x2(i)=eTx(i). Проблема становится

minx,t(igm(i)eTx(i)+t2).               (1)

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

x(i)-x(i-1)-l(i)2k(i)t(i).   (2)

Целевая функция еще не линейна в своих переменных, как требуется для coneprog. Введите новую скалярную переменную y. Заметьте что неравенство t2y эквивалентно неравенству

[2t1-y]1+y                               . (3)

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

minx,t,y(igm(i)eTx(i)+y)                     (4)

подвергните коническим ограничениям на x(i) и t(i) перечисленный в (2) и дополнительное коническое ограничение (3). Коническое ограничение (3) гарантирует это t2y. Поэтому проблема (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 будьте вектором, содержащим все эти переменные.

  • [v(1),v(2)] соответствует 2D переменной x(1).

  • [v(3),v(4)] соответствует 2D переменной x(2).

  • [v(5),v(6)] соответствует 2D переменной x(3).

  • [v(7),v(8)] соответствует 2D переменной x(4).

  • [v(9),v(10)] соответствует 2D переменной x(5).

  • [v(11):v(16)] соответствует 6-D вектору t.

  • v(17) соответствует скалярной переменной y.

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

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

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

x(i)-x(i-1)-l(i)2k(i)t(i).

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

Ascv-bscdscTv-γ.

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

[2t1-y]1+y

путем создания матричного Asc который, когда умножено на v вектор, дает вектор [2t-y]. 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.

Смотрите также

|

Похожие темы

Для просмотра документации необходимо авторизоваться на сайте