Symbolic Computation with Optimisation in MATLAB

Question

So this is a combination symbolic computation and optimisation question. I have a system of three differential equations for phi21, phi31, and phi32. There are four parameters in the equations that I eventually want to optimise for k1, k2, f, and s. I have set up the equations and the construction of a Jacobian in the code below:

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])];

Once I have done this, I am looking to optimise such that f21 = f31 = f32 = 0 and the eigenvalues are all negative. However, I can't figure out how to use my symbolic expressions in some nonlinear optimisation procedure. I have some code that looks like this:

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

where I can specify initial search point, upper bound, and lower bound as well as the vector ceq for my f21, f31, f32 conditions and the vector c for my eigenvalue conditions. There are a couple of problems that I know already. First, the optimisation portion wants variables in the form k(1), k(2), k(3), and k(4) instead of k1, k2, f, and s. Is there a way to easily do this? Second, do I need to convert the symbolic constraints into MATLAB functions? There might be other problems, but I am unsure. Any help would be greatly appreciated :)

Adam Adam · Accepted Answer · 2019-07-05T20:03:34

Use matlabFunction as mentioned by rinkert to convert syms function to function handle
Equality f21 = f31 = f32 is equivalent to f21 - f31 = 0 and f21 - f32 = 0
To get constraints only fulfilled, define the objective function to be a constant function

eq = @(k)0

Please read through the comments

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);

Solution

sol = [0.0116    0.5946   -0.3432    1.0064]

Symbolic Computation with Optimisation in MATLAB

2 Answers