Ограниченная проблема минимизации

Для этой проблемы целевая функция, чтобы минимизировать является простой функцией 2D переменных X и Y:

camxy = @(X,Y)(4 - 2.1.*X.^2 + X.^4./3).*X.^2 + X.*Y + (-4 + 4.*Y.^2).*Y.^2;

Эта функция известна как "бегунок", как описано в L.C.W. Диксон и Г.П. Сзего [1].

Кроме того, проблема имеет нелинейные ограничения и границы.

   x(1)*x(2) + x(1) - x(2) + 1.5 <= 0  (nonlinear constraint)
   10 - x(1)*x(2) <= 0                 (nonlinear constraint)
   0 <= x(1) <= 1                      (bound)
   0 <= x(2) <= 13                     (bound)

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

x1 = linspace(0,1);
y1 = (-x1 - 1.5)./(x1 - 1);
y2 = 10./x1;
[X,Y] = meshgrid(x1,linspace(0,13));
Z = camxy(X,Y);
surf(X,Y,Z,"LineStyle","none")
hold on
z1 = camxy(x1,y1);
z2 = camxy(x1,y2);
plot3(x1,y1,z1,'r-',x1,y2,z2,'r-')
xlim([0 1])
ylim([0 13])
zlim([0,max(Z,[],"all")])
hold off

Создайте переменные оптимизации, проблему и ограничения

Чтобы настроить эту проблему, создайте переменные x оптимизации и y. Установите границы, когда вы создаете переменные.

x = optimvar("x","LowerBound",0,"UpperBound",1);
y = optimvar("y","LowerBound",0,"UpperBound",13);

Создайте цель как выражение оптимизации.

cam = camxy(x,y);

Создайте задачу оптимизации с этой целевой функцией.

prob = optimproblem("Objective",cam);

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

prob.Constraints.cons1 = x*y + x - y + 1.5 <= 0;
prob.Constraints.cons2 = 10 - x*y <= 0;

Рассмотрите проблему.

show(prob)

  OptimizationProblem : 

	Solve for:
       x, y

	minimize :
       (((((4 - (2.1 .* x.^2)) + (x.^4 ./ 3)) .* x.^2) + (x .* y)) + (((-4)
     + (4 .* y.^2)) .* y.^2))


	subject to cons1:
       ((((x .* y) + x) - y) + 1.5) <= 0

	subject to cons2:
       (10 - (x .* y)) <= 0

	variable bounds:
       0 <= x <= 1

       0 <= y <= 13

Установите начальную точку и решите

Установите начальную точку как структуру с полем x равняйтесь 0.5 и y равняйтесь 0.5.

x0.x = 0.5;
x0.y = 0.5;

Решите задачу, задающую patternsearch решатель.

[sol,fval] = solve(prob,x0,"Solver","patternsearch")

Solving problem using patternsearch.
Optimization terminated: mesh size less than options.MeshTolerance
 and constraint violation is less than options.ConstraintTolerance.

sol = struct with fields:
    x: 0.8122
    y: 12.3122

fval = 9.1324e+04

Просмотрите точку решения.

disp(sol.x)

    0.8122

Добавьте визуализацию

Чтобы наблюдать прогресс решателя, задайте опции тот выбор две функции построения графика. Функция построения графика psplotbestf строит лучшее значение целевой функции в каждой итерации и функцию построения графика psplotmaxconstr строит максимальное нарушение ограничений в каждой итерации. Установите эти две функции построения графика в массиве ячеек. Кроме того, информация об отображении о прогрессе решателя в Командном окне путем установки Display опция к 'iter'.

options = optimoptions(@patternsearch,...
    "PlotFcn",{@psplotbestf,@psplotmaxconstr},...
    "Display","iter");

Запустите решатель, включая options аргумент.

[sol,fval] = solve(prob,x0,"Solver","patternsearch","Options",options)

Solving problem using patternsearch.

                                      Max
  Iter   Func-count       f(x)      Constraint   MeshSize      Method
    0         1     0.373958         9.75       0.9086    
    1        18       113581    1.617e-10        0.001   Increase penalty
    2       148        92267            0        1e-05   Increase penalty
    3       374      91333.2            0        1e-07   Increase penalty
    4       639        91324            0        1e-09   Increase penalty
Optimization terminated: mesh size less than options.MeshTolerance
 and constraint violation is less than options.ConstraintTolerance.

sol = struct with fields:
    x: 0.8122
    y: 12.3122

fval = 9.1324e+04

Нелинейные ограничения вызывают patternsearch решить много подпроблем в каждой итерации. Как показано и в графиках и в итеративном отображении, процесс решения имеет немного итераций. Однако Func-count столбец в итеративном отображении показывает много вычислений функции на итерацию. И графики и итеративное отображение показывают, что начальная точка неосуществима, и что целевая функция является низкой в начальной точке. Во время процесса решения значение целевой функции первоначально увеличивается, затем уменьшается к своему окончательному значению.

Неподдерживаемые функции

Если ваши объективные или нелинейные ограничительные функции не являются Поддерживаемыми Операциями для Переменных и выражений Оптимизации, используйте fcn2optimexpr преобразовывать их в форму, подходящую для подхода, основанного на проблеме. Например, предположите это вместо ограничения $$xy\ge 10$$ у вас есть ограничение $$I_1(x)+I_1(y)\ge 10$$, где $$I_1(x)$$ модифицированная Функция Бесселя besseli(1,x). (Функции Бесселя не являются поддерживаемыми функциями.) Создают это ограничение с помощью fcn2optimexpr можно следующим образом. Сначала создайте выражение оптимизации для $$I_1(x)+I_1(y)$$.

bfun = fcn2optimexpr(@(t,u)besseli(1,t) + besseli(1,u),x,y);

Затем замените ограничение cons2 с ограничением bfun >= 10.

prob.Constraints.cons2 = bfun >= 10;

Решите задачу. Решение отличается, потому что область ограничений отличается.

[sol2,fval2] = solve(prob,x0,"Solver","patternsearch","Options",options)

Solving problem using patternsearch.

                                      Max
  Iter   Func-count       f(x)      Constraint   MeshSize      Method
    0         1     0.373958        9.484       0.9307    
    1        18       113581            0        0.001   Increase penalty
    2        78      962.841            0        1e-05   Increase penalty
    3       394      960.942            0        1e-07   Increase penalty
    4       531       960.94            0    8.511e-15   Update multipliers
Optimization terminated: mesh size less than options.MeshTolerance
 and constraint violation is less than options.ConstraintTolerance.

sol2 = struct with fields:
    x: 0.4998
    y: 3.9981

fval2 = 960.9401

Ссылки

[1] Диксон, L. C. W. и G.P. Szego (редакторы).. К глобальной оптимизации 2. Северная Голландия: Elsevier Science Ltd., Амстердам, 1978.