在 MATLAB 中使用三角函数最小化方程和约束

Question

我有一个优化脚本（如下），它给了我错误：

"未定义函数 'cos' 用于输入参数类型 'optim.problemdef.OptimizationExpression'。”

为了确保cos 函数没有问题，我将confn2 中的cos 更改为sin，然后得到：

"未定义函数'sin'用于输入参数类型 'optim.problemdef.OptimizationExpression'。”

但是，如果我在命令行窗口中输入 cos(pi)，我会得到 -1。

我之前已经成功运行过这个脚本，除了 sigma（出现在 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，因为 confn2 中有 cos(sigma)。从 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