Символьные вычисления с оптимизацией в MATLAB - PullRequest
Новая
       33

Символьные вычисления с оптимизацией в MATLAB

1 голос
/ 05 июля 2019

Так что это вопрос сочетания символических вычислений и оптимизации. У меня есть система из трех дифференциальных уравнений для 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

, где я могу указать начальную точку поиска, верхнюю границу и нижнюю границу, а также вектор ceq для моих условий f21, f31, f32 и вектор c для моих условий на собственные значения. Есть пара проблем, которые я уже знаю. Во-первых, для части оптимизации нужны переменные в виде k(1), k(2), k(3) и k(4) вместо k1, k2, f и s. Есть ли способ легко сделать это? Во-вторых, нужно ли преобразовывать символические ограничения в функции MATLAB? Могут быть и другие проблемы, но я не уверен. Любая помощь будет принята с благодарностью:)

Ответы [ 2 ]

1 голос
/ 05 июля 2019
  • Используйте matlabFunction, как указано rinkert , чтобы преобразовать syms функция до function handle
  • Равенство f21 = f31 = f32 эквивалентно f21 - f31 = 0 and f21 - f32 = 0
  • Чтобы только ограничения выполнялись, определите целевую функцию, которая будет постоянная функция 
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]
1 голос
/ 05 июля 2019

Вы преобразуете все необходимые символьные выражения (индивидуально, f21, f31, f32 или seq) в исполняемые функции Matlab, используя matlabFunction. Это сделает их исполняемыми (и выведет двойные вместо символических значений), и позволит им принимать несколько входных аргументов, отсортированных по алфавиту.

Таким образом, matlabFunction(seq) приведет к анонимной функции, которая принимает (f,k1,k2,s) в качестве входных аргументов. Вы также можете сохранить функции в файле, используя аргумент 'File' matlabFunction.

Чтобы позволить этой функции принимать вектор параметров, которые вы хотите оптимизировать, вы можете написать небольшую «обертку»: 

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

Я бы порекомендовал записать сохранить все функции в файлах (также в качестве оболочки), поскольку объектная функция имеет свое собственное рабочее пространство (то есть не может получить доступ к дескрипторам анонимных функций в базовом рабочем пространстве).

...