Ступенчатая консольная проблема проектирования луча

Ступенчатый консольный луч поддерживается в одном конце, и загрузка применяется в свободном конце как показано в следующем рисунке. Луч должен смочь поддержать данную загрузку $$P$$ на фиксированном расстоянии $$L$$ от поддержки. Разработчики луча могут варьироваться ширина ($$b_i$$) и высота ($$h_i$$) из каждого шага или раздела. Каждый раздел консоли имеет ту же длину, $$l=L_1$$.

Объем луча

Объем луча $$V$$ сумма объема отдельных разделов.

$$V=l(b_1{h}_1+b_2{h}_2+b_3{h}_3+b_4{h}_4+b_5{h}_5)$$.

Ограничения на проект: изгиб напряжения

Рассмотрите один консольный луч с центром координат в центре его сечения в свободном конце луча. Изгибающееся напряжение в точке $$(x,y,z)$$ в луче дан уравнением

$${\sigma }_b=M(x)y/I$$,

где $$M(x)$$ изгибающий момент в $$x$$, $$x$$ расстояние от загрузки конца, и $$I$$ момент области инерции луча.

В ступенчатом консольном показанном на рисунке луче максимальный момент каждого раздела луча $$P{D}_i$$, где $$D_i$$ максимальное расстояние от загрузки конца $$P$$ для каждого раздела луча. Поэтому максимальное напряжение для i-ого раздела луча $${\sigma }_i$$ дают

$${\sigma }_i=P{D}_i(h_i/2)/I_i$$,

где максимальное напряжение появляется в ребре луча, $$y=h_i/2$$. Моментом области инерции i-ого раздела луча дают

$$I_i=b_i{h}_i^3/12$$.

Замена этим выражением в уравнение для $${\sigma }_i$$ дает

$${\sigma }_i=6P{D}_i/b_i{h}_i^2$$.

Изгибающееся напряжение в каждой части консоли не должно превышать максимальное допустимое напряжение $${\sigma }_{max}$$. Поэтому пять изгибающихся ограничений напряжения (один для каждого шага консоли):

Ограничения на проект: закончите отклонение

Можно вычислить отклонение конца консоли с помощью второй теоремы Кастиглиано, которая утверждает это

$$\delta =\frac{\partial U}{\partial P}$$,

где $$\delta$$ отклонение луча, и $$U$$ энергия, сохраненная в луче из-за приложенной силы $$P$$.

Энергией, сохраненной в консольном луче, дают

$$U={\int }_0^L{M}^2/2EIdx$$,

где $$M$$ момент приложенной силы в $$x$$.

Учитывая, что $$M=Px$$ для консольного луча можно записать предыдущее уравнение как

$$U=P^2/2E{\int }_0^l[(x+4l{)}^2/I_1+(x+3l{)}^2/I_2+(x+2l{)}^2/I_3+(x+l{)}^2/I_4+x^2/I_5]dx$$,

где $$I_n$$ момент области инерции энной части консоли. Оценка интеграла дает это выражение для $$U$$

$$U=(P^2/2)(l^3/3E)(61/I_1+37/I_2+19/I_3+7/I_4+1/I_5)$$.

Применяя теорему Кастиглиано, отклонением конца луча дают

$$\delta =P{l}^3/3E(61/I_1+37/I_2+19/I_3+7/I_4+1/I_5)$$.

Отклонение конца консоли $$\delta$$ должен быть меньше максимального допустимого отклонения $${\delta }_{max}$$, который дает ограничение

$$P{l}^3/3E(61/I_1+37/I_2+19/I_3+7/I_4+1/I_5)\le {\delta }_{max}$$.

Ограничения на проект: Соотношение сторон

Для каждого шага консоли соотношение сторон не должно превышать максимальное допустимое соотношение сторон $$a_{max}$$. Таким образом,

$$h_i/b_i\le a_{max}$$ для $$i=1,...,5$$.

Ограничения на проект: границы и целочисленные ограничения

Первый шаг луча может быть обработан машинным способом к самому близкому сантиметру только. Таким образом, $$b_1$$ и $$h_1$$ должны быть целые числа. Остающиеся переменные непрерывны. Границы на переменных:

Расчетные параметры для этой проблемы

Для проблемы в этом примере загрузка конца, которую должен поддержать луч, $$P=50000N$$.

Длины луча и максимальное отклонение конца:

Максимальное позволенное напряжение на каждом шаге луча $${\sigma }_{max}=14000N/c{m}^2$$.

Модуль молодежи каждого шага луча $$E=2\times 10^{7}N/c{m}^2$$.

Основанный на проблеме Setup

Создайте переменные оптимизации для этой проблемы. Переменные ширины и высоты для первого раздела луча имеют целое число типа, таким образом, необходимо создать их отдельно из других четырех переменных, которые непрерывны.

b1 = optimvar("b1","Type","integer","LowerBound",1,"UpperBound",5);
h1 = optimvar("h1","Type","integer","LowerBound",30,"UpperBound",65);
bc = optimvar("bc",4,"LowerBound",[2.4 2.4 1 1],"UpperBound",[3.1 3.1 5 5]);
hc = optimvar("hc",4,"LowerBound",[45 45 30 30],"UpperBound",[60 60 65 65]);

Для удобства, помещенного переменные height и width в одну переменные. Можно затем описать цель и ограничения легко в терминах этих переменных.

h = [h1;hc];
b = [b1;bc];

Создайте задачу оптимизации с объемом луча как целевая функция, где каждый шаг (или раздел) луча $$L_1=100$$ cm долго: объем = $$L_1\sum h_i{w}_i$$.

L_1 = 100; % Length of each step of the cantilever
prob = optimproblem("Objective",L_1*dot(h,b));

Создайте ограничения на напряжение.

P = 50000; % End load
E = 2e7; % Young's modulus in N/cm^2
deltaMax = 2.7; % Maximum end deflection
sigmaMax = 14000; % Maximum stress in each section of the beam
aMax = 20; % Maximum aspect ratio in each section of the beam

stress = 6*P*L_1./(b.*(h.^2));
stepnum = (5:-1:1)';
stress = stress.*stepnum;
prob.Constraints.stress = stress <= sigmaMax;

Создайте ограничение на отклонение.

deflectionMultiplier = (P*L_1^3/E)*[244 148 76 28 4];
bh3 = 1./(b.*(h.^3));
prob.Constraints.deflection = deflectionMultiplier*bh3 <= deltaMax;

Создайте ограничения на соотношение сторон.

prob.Constraints.aspect = h <= aMax*b;

Рассмотрите настройку задач.

show(prob)

  OptimizationProblem : 

	Solve for:
       b1, bc, h1, hc
	where:
       b1, h1 integer

	minimize :
       100*h1*b1 + 100*hc(1)*bc(1) + 100*hc(2)*bc(2) + 100*hc(3)*bc(3)
     + 100*hc(4)*bc(4)


	subject to stress:
       arg_LHS <= arg_RHS

       where:

         arg1 = zeros([5, 1]);
         arg1(1) = h1;
         arg1(2:5) = hc;
         arg2 = zeros([5, 1]);
         arg2(1) = b1;
         arg2(2:5) = bc;
         arg_LHS = ((30000000 ./ (arg2(:) .* arg1(:).^2)) .* extraParams{1});
         arg2 = 14000;
         arg1 = arg2(ones(1,5));
         arg_RHS = arg1(:);

         extraParams{1}:

          5
          4
          3
          2
          1



	subject to deflection:
       arg_LHS <= 2.7

       where:

         arg1 = zeros([5, 1]);
         arg1(1) = h1;
         arg1(2:5) = hc;
         arg2 = zeros([5, 1]);
         arg2(1) = b1;
         arg2(2:5) = bc;
         arg_LHS = (extraParams{1} * (1 ./ (arg2(:) .* arg1(:).^3)));

         extraParams{1}:

           610000      370000      190000       70000       10000



	subject to aspect:
       -20*b1 + h1 <= 0
       -20*bc(1) + hc(1) <= 0
       -20*bc(2) + hc(2) <= 0
       -20*bc(3) + hc(3) <= 0
       -20*bc(4) + hc(4) <= 0

	variable bounds:
       1 <= b1 <= 5

       2.4 <= bc(1) <= 3.1
       2.4 <= bc(2) <= 3.1
         1 <= bc(3) <= 5
         1 <= bc(4) <= 5

       30 <= h1 <= 65

       45 <= hc(1) <= 60
       45 <= hc(2) <= 60
       30 <= hc(3) <= 65
       30 <= hc(4) <= 65

Решите задачу

Установите опции использовать умеренную численность населения 150, умеренное максимальное количество поколений 400, немного большее, чем элитное количество по умолчанию 10, маленький функциональный допуск 1e-8 и функция построения графика, показывающая значение функции во время итераций.

opts = optimoptions(@ga, ...
                    'PopulationSize', 150, ...
                    'MaxGenerations', 400, ...
                    'EliteCount', 10, ...
                    'FunctionTolerance', 1e-8, ...
                    'PlotFcn', @gaplotbestf);

Решите задачу, задав ga решатель и опции.

rng default % For reproducibility
[sol,fval,exitflag] = solve(prob,"Solver","ga","Options",opts)

Solving problem using ga.

Optimization terminated: maximum number of generations exceeded.

sol = struct with fields:
    b1: 3.0000
    bc: [4x1 double]
    h1: 60.0000
    hc: [4x1 double]

fval = 6.3558e+04

exitflag = 
    SolverLimitExceeded

Просмотрите переменные решения и их границы.

widths = [sol.b1;sol.bc];
heights = [sol.h1;sol.hc];
widthbounds = [b1.LowerBound b1.UpperBound;
    bc.LowerBound bc.UpperBound];
heightbounds = [h1.LowerBound h1.UpperBound;
    hc.LowerBound hc.UpperBound];
T = table(widths,heights,widthbounds,heightbounds,...
    'VariableNames',["Width" "Height" "Width Bounds" "Height Bounds"])

T=5×4 table
    Width     Height    Width Bounds    Height Bounds
    ______    ______    ____________    _____________

         3        60       1      5       30    65   
    2.8369    56.632     2.4    3.1       45    60   
    2.6483    50.739     2.4    3.1       45    60   
    2.2205    44.293       1      5       30    65   
    1.7652    35.231       1      5       30    65   

Решение не ни в одной из границ. Первые переменные решения являются оцененным целым числом, как задано.

Добавьте дискретные ограничения переменной нецелого числа

Предположим, что инженеры узнают, что вторым и третьим шагам консоли можно было выбрать ширины и высоты от стандартного набора только. Со сложением этого ограничения проблема идентична той, решенной в [1].

Во-первых, формируйте рисунок дополнительных ограничений, чтобы добавить к оптимизации:

Чтобы решить эту задачу, необходимо задать переменные $$wc(1)$$, $$wc(2)$$, $$hc(1)$$, и $$hc(2)$$ как дискретные переменные. Идеально, вы использовали бы $$S(x_j)$$ как дискретное значение, где $$S$$ представляет допустимое множество значений и $$x_j$$ представляет переменную задачи. Однако вы не можете использовать переменную оптимизации в качестве индекса. Можно обойти эту проблему путем вызова fcn2optimexpr.

widthlist = [2.4,2.6,2.8,3.1];
heightlist = [45 50 55 60];
b23 = optimvar("w23",2,"Type","integer",...
    "LowerBound",1,"UpperBound",length(widthlist));
h23 = optimvar("h23",2,"Type","integer",...
    "LowerBound",1,"UpperBound",length(heightlist));
b45 = optimvar("b45",2,"LowerBound",1,"UpperBound",5);
h45 = optimvar("h45",2,"LowerBound",30,"UpperBound",65);
% Preferred syntax is we = [widthlist(b23(1));widthlist(b23(2))];
% However, this syntax is illegal.
% Instead call fcn2optimexpr.
we = fcn2optimexpr(@(x)[widthlist(x(1));widthlist(x(2))],b23);
he = fcn2optimexpr(@(x)[heightlist(x(1));heightlist(x(2))],h23);

Когда вы сделали ранее, создайте выражения b и h представлять переменные.

b = [b1;we;b45];
h = [h1;he;h45];

Остаток от формулировки задачи эквивалентен ранее.

prob2 = optimproblem("Objective",L_1*dot(h,b));

Создайте ограничения на напряжение.

stress = 6*P*L_1./(b.*(h.^2));
stepnum = (5:-1:1)';
stress = stress.*stepnum;
prob2.Constraints.stress = stress <= sigmaMax;

Создайте ограничение на отклонение.

deflectionMultiplier = (P*L_1^3/E)*[244 148 76 28 4];
bh3 = 1./(b.*(h.^3));
prob2.Constraints.deflection = deflectionMultiplier*bh3 <= deltaMax;

Создайте ограничения на соотношение сторон.

prob2.Constraints.aspect = h <= aMax*b;

Рассмотрите настройку задач.

show(prob2)

  OptimizationProblem : 

	Solve for:
       b1, b45, h1, h23, h45, w23
	where:
       b1, h1, h23, w23 integer

	minimize :
       (100 .* (arg1(:).' * arg2(:)))

       where:

         anonymousFunction2 = @(x)[heightlist(x(1));heightlist(x(2))];
         arg1 = zeros([5, 1]);
         arg1(2:3) = anonymousFunction2(h23);
         arg1(1) = h1;
         arg1(4:5) = h45;
         anonymousFunction1 = @(x)[widthlist(x(1));widthlist(x(2))];
         arg2 = zeros([5, 1]);
         arg2(2:3) = anonymousFunction1(w23);
         arg2(1) = b1;
         arg2(4:5) = b45;


	subject to stress:
       arg_LHS <= arg_RHS

       where:

         anonymousFunction2 = @(x)[heightlist(x(1));heightlist(x(2))];
         arg1 = zeros([5, 1]);
         arg1(2:3) = anonymousFunction2(h23);
         arg1(1) = h1;
         arg1(4:5) = h45;
         anonymousFunction1 = @(x)[widthlist(x(1));widthlist(x(2))];
         arg2 = zeros([5, 1]);
         arg2(2:3) = anonymousFunction1(w23);
         arg2(1) = b1;
         arg2(4:5) = b45;
         arg_LHS = ((30000000 ./ (arg2(:) .* arg1(:).^2)) .* extraParams{1});
         arg2 = 14000;
         arg1 = arg2(ones(1,5));
         arg_RHS = arg1(:);

         extraParams{1}:

          5
          4
          3
          2
          1



	subject to deflection:
       arg_LHS <= 2.7

       where:

         anonymousFunction2 = @(x)[heightlist(x(1));heightlist(x(2))];
         arg1 = zeros([5, 1]);
         arg1(2:3) = anonymousFunction2(h23);
         arg1(1) = h1;
         arg1(4:5) = h45;
         anonymousFunction1 = @(x)[widthlist(x(1));widthlist(x(2))];
         arg2 = zeros([5, 1]);
         arg2(2:3) = anonymousFunction1(w23);
         arg2(1) = b1;
         arg2(4:5) = b45;
         arg_LHS = (extraParams{1} * (1 ./ (arg2(:) .* arg1(:).^3)));

         extraParams{1}:

           610000      370000      190000       70000       10000



	subject to aspect:
       arg_LHS <= arg_RHS

       where:

         arg1 = zeros([5, 1]);
         arg1(1) = h1;
         anonymousFunction2 = @(x)[heightlist(x(1));heightlist(x(2))];
         arg1(2:3) = anonymousFunction2(h23);
         arg1(4:5) = h45;
         arg_LHS = arg1(:);
         anonymousFunction1 = @(x)[widthlist(x(1));widthlist(x(2))];
         arg1 = zeros([5, 1]);
         arg1(2:3) = anonymousFunction1(w23);
         arg1(1) = b1;
         arg1(4:5) = b45;
         arg_RHS = (20 .* arg1(:));

	variable bounds:
       1 <= b1 <= 5

       1 <= b45(1) <= 5
       1 <= b45(2) <= 5

       30 <= h1 <= 65

       1 <= h23(1) <= 4
       1 <= h23(2) <= 4

       30 <= h45(1) <= 65
       30 <= h45(2) <= 65

       1 <= w23(1) <= 4
       1 <= w23(2) <= 4

Решите задачу с дискретными переменными ограничениями

Решите задачу, задав ga решатель и опции.

rng default % For reproducibility
[sol2,fval2,exitflag2] = solve(prob2,"Solver","ga","Options",opts)

Solving problem using ga.

Optimization terminated: maximum number of generations exceeded.

sol2 = struct with fields:
     b1: 3
    b45: [2x1 double]
     h1: 60
    h23: [2x1 double]
    h45: [2x1 double]
    w23: [2x1 double]

fval2 = 6.4803e+04

exitflag2 = 
    SolverLimitExceeded

Объективные повышения стоимости, потому что добавление ограничений может только увеличить цель.

Просмотрите решение и сравните его с его границами.

widths = [sol2.b1;widthlist(sol2.w23(1));widthlist(sol2.w23(2));sol2.b45];
heights = [sol2.h1;heightlist(sol2.h23(1));heightlist(sol2.h23(2));sol2.h45];
widthbounds = [b1.LowerBound b1.UpperBound;
    widthlist(1) widthlist(end);
     widthlist(1) widthlist(end);
    b45.LowerBound b45.UpperBound];
heightbounds = [h1.LowerBound h1.UpperBound;
     heightlist(1) heightlist(end);
     heightlist(1) heightlist(end);
    h45.LowerBound h45.UpperBound];
T = table(widths,heights,widthbounds,heightbounds,...
    'VariableNames',["Width" "Height" "Width Bounds" "Height Bounds"])

T=5×4 table
    Width     Height    Width Bounds    Height Bounds
    ______    ______    ____________    _____________

         3        60       1      5       30    65   
       3.1        55     2.4    3.1       45    60   
       2.6        50     2.4    3.1       45    60   
     2.286     45.72       1      5       30    65   
    1.8532    34.004       1      5       30    65   

Единственная переменная решения, которая является в связанном, является шириной второго раздела, который является 3.1, его максимум.

Ссылки

[1] Thanedar, P. B. и Г. Н. Вэндерплээтс. "Обзор Дискретной Переменной Оптимизации для Проектирования конструкций". Журнал Структурной Разработки 121 (3), 1995, стр 301–306.