Минимизация уравнений и ограничений с помощью тригонометрических функций в MATLAB

Question

У меня есть скрипт оптимизации (ниже), который выдает мне ошибку:

"Неопределенная функция 'cos' для входных аргументов типа 'Optim.problemdef.OptimizationExpression'. "

Просто чтобы убедиться, что это не проблема с функцией cos, я изменил cos в confn2 на sin и затем получил:

"Неопределенная функция 'sin' для входных аргументов типа 'Optim.problemdef.OptimizationExpression'. "

Однако, если я добавлю cos(pi) в командном окне, я получу -1.

Ранее я успешно запускал этот скрипт, за исключением того, что сигма (которая появляется в cos()) не была переменной оптимизации, и, таким образом, функции триггера оценивались, оставляя что-то линейное. Разве этот тип сценария оптимизации не подходит для ограничений с тригонометрическими функциями? Есть ли альтернативы, которые есть в MATLAB?

k1 = optimvar('k1', 'LowerBound', -3, 'UpperBound', 3);
k2 = optimvar('k2', 'LowerBound', -3, 'UpperBound', 3);
f = optimvar('f', 'LowerBound', -3, 'UpperBound', 3);
sigma = optimvar('sigma', 'LowerBound', 0, 'UpperBound', 6.28318530718);

obj = fcn2optimexpr(@eq1, k1, k2, f, sigma);

confn1 = obj == 0;
confn2 = -0.9313129901097519*k1 - 1.4693755421886672*k2 - 0.18532000686683578*f*cos((1/2)*(-0.372795 + 2*sigma)) + 0.9826782255931369*f*sin((1/2)*(-0.372795 + 2*sigma)) <= 0;

prob = optimproblem('Objective', obj);

prob.Constraints.confn1 = confn1;
prob.Constraints.confn2 = confn2;

k0.k1 = 0;
k0.k2 = 0;
k0.f = 0;
k0.sigma = 0;

[sol, fval, exitflag, output] = solve(prob, k0)

function f1 = eq1(k1, k2, f, sigma)
f1 = 0.01308996938995749 - 0.3642198710296203*k1 - 0.6784053942919677*k2 + 0.37064001373367156*f*sin((1/2)*(-0.372795 + 2*sigma));
end

Answer 1

Вместо этого вы можете использовать fmincon, в котором вы также можете указать цель и ограничения.Таким образом, функция cos будет оцениваться численно на каждой итерации и не будет вызываться с типом optim.problemdef.OptimizationExpression.

Вам необходимо указать переменные в виде вектора (начальное предположение x0, нижний предел lb, верхний предел ub), см. Комментарий в коде.

% x0 = [k1 k2 f sigma];
x0 = [0 0 0 0];
lb = [-3 -3 -3 0];
ub = [3 3 3 2*pi];

% options
options = optimoptions('fmincon','Display','iter', 'ConstraintTolerance', 1e-12);

% optimization
x = fmincon(@objfun,x0,[],[],[],[],lb,ub,@constrfun, options)

% check constraints to check fulfilment of constraints
[c, ceq] = constrfun(x)

ВашФункция объекта должна принимать один вектор x, который можно «распаковать» для отдельных переменных.

function f1 = objfun(x)
    k1 = x(1);
    k2 = x(2); 
    f = x(3);
    sigma = x(4);
    f1 = 0.01308996938995749 - 0.3642198710296203*k1 - 0.6784053942919677*k2 + 0.37064001373367156*f*sin((1/2)*(-0.372795 + 2*sigma));
end

Ваши ограничения должны быть определены в функции ограничения, которая возвращает c и ceq, в котором c - это ограничение неравенства c(x) <= 0, а ceq - ограничение равенства ceq(x) = 0.

function [c,ceq] = constrfun(x)
    k1 = x(1);
    k2 = x(2); 
    f = x(3);
    sigma = x(4);
    c = -0.9313129901097519*k1 - 1.4693755421886672*k2 - 0.18532000686683578*f*cos((1/2)*(-0.372795 + 2*sigma)) + 0.9826782255931369*f*sin((1/2)*(-0.372795 + 2*sigma));
    ceq = objfun(x);
end

Answer 2

Вам нужно использовать fcn2optimexpr из-за cos (sigma) в confn2. Начиная с выпуска R2019a, для выражений optimvars разрешены только отношения полиномов.

k1 = optimvar('k1', 'LowerBound', -3, 'UpperBound', 3);
k2 = optimvar('k2', 'LowerBound', -3, 'UpperBound', 3);
f = optimvar('f', 'LowerBound', -3, 'UpperBound', 3);
sigma = optimvar('sigma', 'LowerBound', 0, 'UpperBound', 6.28318530718);

obj = fcn2optimexpr(@eq1, k1, k2, f, sigma);
c2 = fcn2optimexpr(@conexpr2, k1, k2, f, sigma);

confn1 = obj == 0;
confn2 = c2 <= 0;

prob = optimproblem('Objective', obj);

prob.Constraints.confn1 = confn1;
prob.Constraints.confn2 = confn2;

k0.k1 = 0;
k0.k2 = 0;
k0.f = 0;
k0.sigma = 0;

options = optimoptions(prob,'Display','iter', 'ConstraintTolerance', 1e-12);
[sol, fval, exitflag, output] = solve(prob, k0, 'Options',options)

objval = evaluate(obj,sol)
c2val = evaluate(c2,sol)


function f1 = eq1(k1, k2, f, sigma)
f1 = 0.01308996938995749 - 0.3642198710296203*k1 - 0.6784053942919677*k2 + 0.37064001373367156*f*sin((1/2)*(-0.372795 + 2*sigma));
end

function c2 = conexpr2(k1, k2, f, sigma)
c2 = -0.9313129901097519*k1 - 1.4693755421886672*k2 - 0.18532000686683578*f*cos((1/2)*(-0.372795 + 2*sigma)) + 0.9826782255931369*f*sin((1/2)*(-0.372795 + 2*sigma));
end

Как заметил @rinkert, вы используете цель дважды. Нет необходимости иметь объективную функцию. Вы можете столкнуться с проблемой оптимизации только с ограничениями.

k1 = optimvar('k1', 'LowerBound', -3, 'UpperBound', 3);
k2 = optimvar('k2', 'LowerBound', -3, 'UpperBound', 3);
f = optimvar('f', 'LowerBound', -3, 'UpperBound', 3);
sigma = optimvar('sigma', 'LowerBound', 0, 'UpperBound', 6.28318530718);

c1 = fcn2optimexpr(@eq1, k1, k2, f, sigma);
c2 = fcn2optimexpr(@conexpr2, k1, k2, f, sigma);

confn1 = c1 == 0;
confn2 = c2 <= 0;

prob = optimproblem();

prob.Constraints.confn1 = confn1;
prob.Constraints.confn2 = confn2;

k0.k1 = 0;
k0.k2 = 0;
k0.f = 0;
k0.sigma = 0;

options = optimoptions(prob,'Display','iter', 'ConstraintTolerance', 1e-12);
[sol, fval, exitflag, output] = solve(prob, k0, 'Options',options)

c1val = evaluate(c1,sol)
c2val = evaluate(c2,sol)


function f1 = eq1(k1, k2, f, sigma)
f1 = 0.01308996938995749 - 0.3642198710296203*k1 - 0.6784053942919677*k2 + 0.37064001373367156*f*sin((1/2)*(-0.372795 + 2*sigma));
end

function c2 = conexpr2(k1, k2, f, sigma)
c2 = -0.9313129901097519*k1 - 1.4693755421886672*k2 - 0.18532000686683578*f*cos((1/2)*(-0.372795 + 2*sigma)) + 0.9826782255931369*f*sin((1/2)*(-0.372795 + 2*sigma));
end