MATLAB 中的优化符号计算
Symbolic Computation with Optimisation in MATLAB
所以这是一个结合符号计算和优化的问题。我有一个包含 phi21
、phi31
和 phi32
的三个微分方程组。方程式中有四个参数我最终想优化 k1
、k2
、f
和 s
。我在下面的代码中设置了方程式和雅可比行列式的构造:
syms phi21 phi31 phi32 k1 k2 f s a
w1 = (2*pi/24)*0.99;
w2 = (2*pi/24)*1.01;
w3 = (2*pi/24)*1.02;
f21 = (w2 - w1) + (1/3)*(-2*k1*sin(phi21) - 2*k2*sin(2*phi21) - k1*sin(phi31) - k2*sin(2*phi31) + k1*sin(phi32) + k2*sin(2*phi32)) + 2*f*cos((2*s - phi21)/2)*sin(-phi21/2);
f31 = (w3 - w1) + (1/3)*(k1*sin(phi21) + k2*sin(2*phi21) - k1*sin(phi31) - k2*sin(2*phi31) - 2*k1*sin(phi32) - 2*k2*sin(2*phi32)) + 2*f*cos((2*s - phi31)/2)*sin(-phi31/2);
f32 = (w3 - w2) + (1/3)*(-k1*sin(phi21) - k2*sin(2*phi21) - 2*k1*sin(phi31) - 2*k2*sin(2*phi31) - k1*sin(phi32) - k2*sin(2*phi32)) + + 2*f*cos((2*s - phi32)/2)*sin(-phi32/2);
df21d21 = diff(f21, phi21);
df21d31 = diff(f21, phi31);
df21d32 = diff(f21, phi32);
df31d21 = diff(f31, phi21);
df31d31 = diff(f31, phi31);
df31d32 = diff(f31, phi32);
df32d21 = diff(f32, phi21);
df32d31 = diff(f32, phi31);
df32d32 = diff(f32, phi32);
J = [df21d21 df21d31 df21d32; df31d21 df31d31 df31d32; df32d21 df32d31 df32d32];
lambda = eig(J);
rlambda = real(lambda);
srlambda = subs(rlambda, [phi21, phi31, phi32], [0.35475, 0.58305, 0.2271]);
seq = [subs(f21, [phi21, phi31, phi32], [0.35475, 0.58305, 0.2271]), subs(f31, [phi21, phi31, phi32], [0.35475, 0.58305, 0.2271]), subs(f32, [phi21, phi31, phi32], [0.35475, 0.58305, 0.2271])];
完成此操作后,我希望进行优化以使 f21 = f31 = f32 = 0 并且特征值均为负数。但是,我不知道如何在某些非线性优化过程中使用我的符号表达式。我有一些代码如下所示:
x0 = [];
lb = [];
ub = [];
[sol, fval, exitflag, output] = fmincon(@eq1, x0, A, b, Aeq, beq, lb, ub, @constraints)
function objfun = eq1(k)
objfun = ;
end
function [c, ceq] = constraints(k)
c = [];
ceq = [];
end
我可以在其中为我的 f21
、f31
、f32
条件指定初始搜索点、上限和下限以及向量 ceq
和我的特征值条件的向量 c
。
我已经知道有几个问题。首先,优化部分需要 k(1)
、k(2)
、k(3)
和 k(4)
形式的变量,而不是 k1
、k2
、f
,和 s
。有没有办法轻松做到这一点?其次,我是否需要将符号约束转换为 MATLAB 函数?可能还有其他问题,但我不确定。任何帮助将不胜感激:)
您使用 matlabFunction
将所有需要的符号表达式(分别为 f21
、f31
、f32
或 seq
)转换为可执行的 Matlab 函数.这将使它们可执行(并输出双精度而不是符号值),并让它们接受按字母顺序排序的多个输入参数。
因此,matlabFunction(seq)
将生成一个以 (f,k1,k2,s)
作为输入参数的匿名函数。您还可以使用 matlabFunction
的 'File'
参数将函数存储在文件中。
为了让这个函数接受一个你想要优化的参数向量,你可以写一个小的 'wrapper':
f_seq = matlabFunction(seq); % anonymous function that takes (f,k1,k2,s) as inputs
f_seq2 = @(p) f_seq(p(1),p(2),p(3),p(4)); % anonymous function that takes a 4 element vector
我建议将所有函数保存在文件中(也是包装器),因为对象函数有自己的工作区(即无法访问基础工作区中的匿名函数句柄)。
- 使用rinkert提到的
matlabFunction
转换syms
功能 function handle
- 等于
f21 = f31 = f32
等价于f21 - f31 = 0 and f21 - f32 = 0
- 要仅满足约束条件,请将 objective 函数定义为
一个常数函数
eq = @(k)0
请仔细阅读评论
syms phi21 phi31 phi32 k1 k2 f s a
w1 = (2*pi/24)*0.99;
w2 = (2*pi/24)*1.01;
w3 = (2*pi/24)*1.02;
f21 = (w2 - w1) + (1/3)*(-2*k1*sin(phi21) - 2*k2*sin(2*phi21) - k1*sin(phi31) - k2*sin(2*phi31) + k1*sin(phi32) + k2*sin(2*phi32)) + 2*f*cos((2*s - phi21)/2)*sin(-phi21/2);
f31 = (w3 - w1) + (1/3)*(k1*sin(phi21) + k2*sin(2*phi21) - k1*sin(phi31) - k2*sin(2*phi31) - 2*k1*sin(phi32) - 2*k2*sin(2*phi32)) + 2*f*cos((2*s - phi31)/2)*sin(-phi31/2);
f32 = (w3 - w2) + (1/3)*(-k1*sin(phi21) - k2*sin(2*phi21) - 2*k1*sin(phi31) - 2*k2*sin(2*phi31) - k1*sin(phi32) - k2*sin(2*phi32)) + + 2*f*cos((2*s - phi32)/2)*sin(-phi32/2);
df21d21 = diff(f21, phi21);
df21d31 = diff(f21, phi31);
df21d32 = diff(f21, phi32);
df31d21 = diff(f31, phi21);
df31d31 = diff(f31, phi31);
df31d32 = diff(f31, phi32);
df32d21 = diff(f32, phi21);
df32d31 = diff(f32, phi31);
df32d32 = diff(f32, phi32);
J = [df21d21 df21d31 df21d32; df31d21 df31d31 df31d32; df32d21 df32d31 df32d32];
lambda = eig(J);
rlambda = real(lambda);
srlambda = subs(rlambda, [phi21, phi31, phi32], [0.35475, 0.58305, 0.2271]);
seq = [subs(f21, [phi21, phi31, phi32], [0.35475, 0.58305, 0.2271]), subs(f31, [phi21, phi31, phi32], [0.35475, 0.58305, 0.2271]), subs(f32, [phi21, phi31, phi32], [0.35475, 0.58305, 0.2271])];
% Transform syms function to function handles
f21 = matlabFunction(seq(1));
f31 = matlabFunction(seq(2));
f32 = matlabFunction(seq(3));
lambda = matlabFunction(srlambda);
% Inequality constraint, input is passed as a vector
c = @(k)lambda(k(1), k(2), k(3), k(4));
% Equality constraint, input is passed as a vector
% f21 = f31 = f32 --> f21 -f31 = 0 and f21 -f32 = 0
ceq = @(k)[f21(k(1), k(2), k(3), k(4))-f31(k(1), k(2), k(3), k(4));...
f21(k(1), k(2), k(3), k(4))-f32(k(1), k(2), k(3), k(4))];
% Combine all the constraints to one function handle
constraints = @(k)deal(c(k),ceq(k));
% Only need the constraints to be satisfied, define a constant objective
% function
eq1 = @(k)0;
% A random starting guess, lower bound, upper bound
% You can change this part to what you want
x0 = ones(1,4);
lb = [-inf, -inf, -inf, -inf];
ub = [inf, inf, inf, inf];
% No linear constraints
A = [];
b = [];
Aeq = [];
beq = [];
[sol, fval, exitflag, output] = fmincon(eq1, x0, A, b, Aeq, beq, lb, ub, constraints);
解决方案
sol = [0.0116 0.5946 -0.3432 1.0064]
所以这是一个结合符号计算和优化的问题。我有一个包含 phi21
、phi31
和 phi32
的三个微分方程组。方程式中有四个参数我最终想优化 k1
、k2
、f
和 s
。我在下面的代码中设置了方程式和雅可比行列式的构造:
syms phi21 phi31 phi32 k1 k2 f s a
w1 = (2*pi/24)*0.99;
w2 = (2*pi/24)*1.01;
w3 = (2*pi/24)*1.02;
f21 = (w2 - w1) + (1/3)*(-2*k1*sin(phi21) - 2*k2*sin(2*phi21) - k1*sin(phi31) - k2*sin(2*phi31) + k1*sin(phi32) + k2*sin(2*phi32)) + 2*f*cos((2*s - phi21)/2)*sin(-phi21/2);
f31 = (w3 - w1) + (1/3)*(k1*sin(phi21) + k2*sin(2*phi21) - k1*sin(phi31) - k2*sin(2*phi31) - 2*k1*sin(phi32) - 2*k2*sin(2*phi32)) + 2*f*cos((2*s - phi31)/2)*sin(-phi31/2);
f32 = (w3 - w2) + (1/3)*(-k1*sin(phi21) - k2*sin(2*phi21) - 2*k1*sin(phi31) - 2*k2*sin(2*phi31) - k1*sin(phi32) - k2*sin(2*phi32)) + + 2*f*cos((2*s - phi32)/2)*sin(-phi32/2);
df21d21 = diff(f21, phi21);
df21d31 = diff(f21, phi31);
df21d32 = diff(f21, phi32);
df31d21 = diff(f31, phi21);
df31d31 = diff(f31, phi31);
df31d32 = diff(f31, phi32);
df32d21 = diff(f32, phi21);
df32d31 = diff(f32, phi31);
df32d32 = diff(f32, phi32);
J = [df21d21 df21d31 df21d32; df31d21 df31d31 df31d32; df32d21 df32d31 df32d32];
lambda = eig(J);
rlambda = real(lambda);
srlambda = subs(rlambda, [phi21, phi31, phi32], [0.35475, 0.58305, 0.2271]);
seq = [subs(f21, [phi21, phi31, phi32], [0.35475, 0.58305, 0.2271]), subs(f31, [phi21, phi31, phi32], [0.35475, 0.58305, 0.2271]), subs(f32, [phi21, phi31, phi32], [0.35475, 0.58305, 0.2271])];
完成此操作后,我希望进行优化以使 f21 = f31 = f32 = 0 并且特征值均为负数。但是,我不知道如何在某些非线性优化过程中使用我的符号表达式。我有一些代码如下所示:
x0 = [];
lb = [];
ub = [];
[sol, fval, exitflag, output] = fmincon(@eq1, x0, A, b, Aeq, beq, lb, ub, @constraints)
function objfun = eq1(k)
objfun = ;
end
function [c, ceq] = constraints(k)
c = [];
ceq = [];
end
我可以在其中为我的 f21
、f31
、f32
条件指定初始搜索点、上限和下限以及向量 ceq
和我的特征值条件的向量 c
。
我已经知道有几个问题。首先,优化部分需要 k(1)
、k(2)
、k(3)
和 k(4)
形式的变量,而不是 k1
、k2
、f
,和 s
。有没有办法轻松做到这一点?其次,我是否需要将符号约束转换为 MATLAB 函数?可能还有其他问题,但我不确定。任何帮助将不胜感激:)
您使用 matlabFunction
将所有需要的符号表达式(分别为 f21
、f31
、f32
或 seq
)转换为可执行的 Matlab 函数.这将使它们可执行(并输出双精度而不是符号值),并让它们接受按字母顺序排序的多个输入参数。
因此,matlabFunction(seq)
将生成一个以 (f,k1,k2,s)
作为输入参数的匿名函数。您还可以使用 matlabFunction
的 'File'
参数将函数存储在文件中。
为了让这个函数接受一个你想要优化的参数向量,你可以写一个小的 'wrapper':
f_seq = matlabFunction(seq); % anonymous function that takes (f,k1,k2,s) as inputs
f_seq2 = @(p) f_seq(p(1),p(2),p(3),p(4)); % anonymous function that takes a 4 element vector
我建议将所有函数保存在文件中(也是包装器),因为对象函数有自己的工作区(即无法访问基础工作区中的匿名函数句柄)。
- 使用rinkert提到的
matlabFunction
转换syms
功能function handle
- 等于
f21 = f31 = f32
等价于f21 - f31 = 0 and f21 - f32 = 0
- 要仅满足约束条件,请将 objective 函数定义为 一个常数函数
eq = @(k)0
请仔细阅读评论
syms phi21 phi31 phi32 k1 k2 f s a
w1 = (2*pi/24)*0.99;
w2 = (2*pi/24)*1.01;
w3 = (2*pi/24)*1.02;
f21 = (w2 - w1) + (1/3)*(-2*k1*sin(phi21) - 2*k2*sin(2*phi21) - k1*sin(phi31) - k2*sin(2*phi31) + k1*sin(phi32) + k2*sin(2*phi32)) + 2*f*cos((2*s - phi21)/2)*sin(-phi21/2);
f31 = (w3 - w1) + (1/3)*(k1*sin(phi21) + k2*sin(2*phi21) - k1*sin(phi31) - k2*sin(2*phi31) - 2*k1*sin(phi32) - 2*k2*sin(2*phi32)) + 2*f*cos((2*s - phi31)/2)*sin(-phi31/2);
f32 = (w3 - w2) + (1/3)*(-k1*sin(phi21) - k2*sin(2*phi21) - 2*k1*sin(phi31) - 2*k2*sin(2*phi31) - k1*sin(phi32) - k2*sin(2*phi32)) + + 2*f*cos((2*s - phi32)/2)*sin(-phi32/2);
df21d21 = diff(f21, phi21);
df21d31 = diff(f21, phi31);
df21d32 = diff(f21, phi32);
df31d21 = diff(f31, phi21);
df31d31 = diff(f31, phi31);
df31d32 = diff(f31, phi32);
df32d21 = diff(f32, phi21);
df32d31 = diff(f32, phi31);
df32d32 = diff(f32, phi32);
J = [df21d21 df21d31 df21d32; df31d21 df31d31 df31d32; df32d21 df32d31 df32d32];
lambda = eig(J);
rlambda = real(lambda);
srlambda = subs(rlambda, [phi21, phi31, phi32], [0.35475, 0.58305, 0.2271]);
seq = [subs(f21, [phi21, phi31, phi32], [0.35475, 0.58305, 0.2271]), subs(f31, [phi21, phi31, phi32], [0.35475, 0.58305, 0.2271]), subs(f32, [phi21, phi31, phi32], [0.35475, 0.58305, 0.2271])];
% Transform syms function to function handles
f21 = matlabFunction(seq(1));
f31 = matlabFunction(seq(2));
f32 = matlabFunction(seq(3));
lambda = matlabFunction(srlambda);
% Inequality constraint, input is passed as a vector
c = @(k)lambda(k(1), k(2), k(3), k(4));
% Equality constraint, input is passed as a vector
% f21 = f31 = f32 --> f21 -f31 = 0 and f21 -f32 = 0
ceq = @(k)[f21(k(1), k(2), k(3), k(4))-f31(k(1), k(2), k(3), k(4));...
f21(k(1), k(2), k(3), k(4))-f32(k(1), k(2), k(3), k(4))];
% Combine all the constraints to one function handle
constraints = @(k)deal(c(k),ceq(k));
% Only need the constraints to be satisfied, define a constant objective
% function
eq1 = @(k)0;
% A random starting guess, lower bound, upper bound
% You can change this part to what you want
x0 = ones(1,4);
lb = [-inf, -inf, -inf, -inf];
ub = [inf, inf, inf, inf];
% No linear constraints
A = [];
b = [];
Aeq = [];
beq = [];
[sol, fval, exitflag, output] = fmincon(eq1, x0, A, b, Aeq, beq, lb, ub, constraints);
解决方案
sol = [0.0116 0.5946 -0.3432 1.0064]