Matrix Optimization using optimization toolbox - "Objective must be a scalar OptimizationExpression or a struct containing a scalar OptimizationExpression."

38 views (last 30 days)
Samuele Bolotta on 13 Feb 2021
Commented: Samuele Bolotta on 13 Feb 2021
Hello everyone,
My task is to find x and y so that if I multiply them by G_max_chl2 and G_max_glu2 (two scalar values that I pre-defined in previous parts of the code) respectively, this:
((G_max_chl2 * x) .* (1 - exp(-t / tau_rise_In2)) .* (exp(-t / tau_decay_In2)) * (Vm - EChl2)) + ((G_max_glu2 * y) .* (1 - exp(-t / tau_rise_Ex2)) .* exp(-t / tau_decay_Ex2) * (Vm - EGlu2))
becomes the same as this:
((G_max_chl) .* (1 - exp(-t / tau_rise_In)) .* (exp(-t / tau_decay_In)) * (Vm - EChl)) + (G_max_glu .* (1 - exp(-t / tau_rise_Ex)) .* exp(-t / tau_decay_Ex) * (Vm - EGlu))
1 All the variables above are the same for the two equations, except for the four "G_max" values. The only two variables that change are G_max_chl2 and G_max_glu2, indeed.
2 Some of these variables are matrices, but again they are not affected by the optimization problem and should remain the same, while the two scalars G_max_chl2 and G_max_glu2 should change
prob = optimproblem('ObjectiveSense','min');
x = optimvar('x',1);
y = optimvar('y',1);
expr = ((G_max_chl2 * x) .* (1 - exp(-t / tau_rise_In2)) .* ...
(exp(-t / tau_decay_In2)) * (Vm - EChl2)) + ((G_max_glu2 * y) .* ...
(1 - exp(-t / tau_rise_Ex2)) .* exp(-t / tau_decay_Ex2) * (Vm - EGlu2)) == ((G_max_chl) .* ...
(1 - exp(-t / tau_rise_In)) .* (exp(-t / tau_decay_In)) * ...
(Vm - EChl)) + (G_max_glu .* (1 - exp(-t / tau_rise_Ex)) .* exp(-t / tau_decay_Ex) * (Vm - EGlu));
prob.Objective = expr;
% Create constraints in the problem
cons1 = ((G_max_chl2 * x) .* (1 - exp(-t / tau_rise_In2)) .* ...
(exp(-t / tau_decay_In2)) * (Vm - EChl2)) + ((G_max_glu2 * y) .* ...
(1 - exp(-t / tau_rise_Ex2)) .* exp(-t / tau_decay_Ex2) * (Vm - EGlu2)) - ((G_max_chl) .* ...
(1 - exp(-t / tau_rise_In)) .* (exp(-t / tau_decay_In)) * ...
(Vm - EChl)) + (G_max_glu .* (1 - exp(-t / tau_rise_Ex)) .* exp(-t / tau_decay_Ex) * (Vm - EGlu)) < 0.1;
cons2 = ((G_max_chl2 * x) .* (1 - exp(-t / tau_rise_In2)) .* ...
(exp(-t / tau_decay_In2)) * (Vm - EChl2)) + ((G_max_glu2 * y) .* ...
(1 - exp(-t / tau_rise_Ex2)) .* exp(-t / tau_decay_Ex2) * (Vm - EGlu2)) - ((G_max_chl) .* ...
(1 - exp(-t / tau_rise_In)) .* (exp(-t / tau_decay_In)) * ...
(Vm - EChl)) + (G_max_glu .* (1 - exp(-t / tau_rise_Ex)) .* exp(-t / tau_decay_Ex) * (Vm - EGlu)) > -0.1;
prob.Constraints.cons1 = cons1;
prob.Constraints.cons2 = cons2;
[sol,fval,exitflag,output] = solve(prob);
I am not sure what is not working. Thanks!

Matt J on 13 Feb 2021
Edited: Matt J on 13 Feb 2021
You do not have a minimization problem. You just have a system of equations that you are trying to solve. For that, you would use the EquationProblem framework.
x = optimvar('x',1);
y = optimvar('y',1);
expr = ((G_max_chl2 * x) .* (1 - exp(-t / tau_rise_In2)) .* ...
(exp(-t / tau_decay_In2)) * (Vm - EChl2)) + ((G_max_glu2 * y) .* ...
(1 - exp(-t / tau_rise_Ex2)) .* exp(-t / tau_decay_Ex2) * (Vm - EGlu2)) == ((G_max_chl) .* ...
(1 - exp(-t / tau_rise_In)) .* (exp(-t / tau_decay_In)) * ...
(Vm - EChl)) + (G_max_glu .* (1 - exp(-t / tau_rise_Ex)) .* exp(-t / tau_decay_Ex) * (Vm - EGlu));
% Create constraints in the problem
cons1 = ((G_max_chl2 * x) .* (1 - exp(-t / tau_rise_In2)) .* ...
(exp(-t / tau_decay_In2)) * (Vm - EChl2)) + ((G_max_glu2 * y) .* ...
(1 - exp(-t / tau_rise_Ex2)) .* exp(-t / tau_decay_Ex2) * (Vm - EGlu2)) - ((G_max_chl) .* ...
(1 - exp(-t / tau_rise_In)) .* (exp(-t / tau_decay_In)) * ...
(Vm - EChl)) + (G_max_glu .* (1 - exp(-t / tau_rise_Ex)) .* exp(-t / tau_decay_Ex) * (Vm - EGlu)) < 0.1;
cons2 = ((G_max_chl2 * x) .* (1 - exp(-t / tau_rise_In2)) .* ...
(exp(-t / tau_decay_In2)) * (Vm - EChl2)) + ((G_max_glu2 * y) .* ...
(1 - exp(-t / tau_rise_Ex2)) .* exp(-t / tau_decay_Ex2) * (Vm - EGlu2)) - ((G_max_chl) .* ...
(1 - exp(-t / tau_rise_In)) .* (exp(-t / tau_decay_In)) * ...
(Vm - EChl)) + (G_max_glu .* (1 - exp(-t / tau_rise_Ex)) .* exp(-t / tau_decay_Ex) * (Vm - EGlu)) > -0.1;
prob = eqnproblem;
prob.Equations.eqn1=expr;
prob.Equations.eqn2 = cons1;
prob.Equations.eqn3 = cons2;
sol= solve(prob);
Samuele Bolotta on 13 Feb 2021
Yes, this is what I did:
% EQUATION
% Create a nonlinear equation
x = optimvar('x',1);
y = optimvar('y',1);
% Define function to minimize
eq1 = ((G_max_chl2 * x) .* ((1 - exp(-t / tau_rise_In2)) .* ...
(exp(-t / tau_decay_In2))) * (Vm - EChl2)) + ((G_max_glu2 * y .* ...
((1 - exp(-t / tau_rise_Ex2)) .* exp(-t / tau_decay_Ex2))* ...
(Vm - EGlu2))) == ((G_max_chl) .* ((1 - exp(-t / tau_rise_In)) .* ...
exp(-t / tau_decay_In)) * (Vm - EChl)) + ((G_max_glu) .* ...
((1 - exp(-t / tau_rise_Ex)) .* exp(-t / tau_decay_Ex)) * (Vm - EGlu));
% Create an equation problem, and place the equation in the problem
prob = eqnproblem;
prob.Equations.eq1 = eq1;
% Show the problem
show(prob);
% Specify the initial point as a structure
x0.x = G_max_chl / G_max_chl2;
x0.y = G_max_glu / G_max_glu2;
[sol,fval,exitflag] = solve(prob,x0);
% View the solution point and convert to double
disp(sol.x);
disp(sol.y);
x = sol.x;
y = sol.y;
% Update maximal conductances
G_max_chl2_new = G_max_chl2 * x;
G_max_glu2_new = G_max_glu2 * y;
% Update conductances CPSC2
Gi2_new = G_max_chl2_new .* gat_I2; % Inhibitory
Ge2_new = G_max_glu2_new .* gat_E2; % Excitatory
% Update CPSC2
IPSC2_new = Gi2_new * (Vm - EChl2); %Inhibitory
EPSC2_new = Ge2_new * (Vm - EGlu2); %Excitatory
CPSC2_new = IPSC2_new + EPSC2_new; %Compound
% Update difference between currents now
New_Diff_CPSC = abs((((G_max_chl2_new) .* ((1 - exp(-t / tau_rise_In2)) .* ...
(exp(-t / tau_decay_In2))) * (Vm - EChl2)) + ((G_max_glu2_new) .* ...
((1 - exp(-t / tau_rise_Ex2)) .* exp(-t / tau_decay_Ex2)) * (Vm - EGlu2))) - (((G_max_chl) .* ...
((1 - exp(-t / tau_rise_In)) .* exp(-t / tau_decay_In)) * ...
(Vm - EChl)) + ((G_max_glu) .* ((1 - exp(-t / tau_rise_Ex)) .* exp(-t / tau_decay_Ex)) * (Vm - EGlu))));
Thank you again,
Sam