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

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!

Answer 1

[Matt J]

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]

Hi Matt, thank you very much for your help!

When I try to use the following code (what you suggested with some parentheses modified because I noticed I had put some in the wrong places):

% OPTIMIZATION PROBLEM
% Create a linear programming problem for minimization

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

The output is:

Error using optimization (line 109)
Equation must be an OptimizationEquality or a struct containing OptimizationEquality.

Line 109 is this one:

prob.Equations.eqn1=expr;

I think the problem is I have matrices in that code, and so I should put them in a struct maybe? Just to give you an idea, if you want to run this yourself, this is all the data of the simulation:

%Set CPSC1
G_max_chl = 20;          % Maximal conductance of chloride
G_max_glu = 30;          % Maximal conductance of glutamate
Vm = -50;           % Array of holding potentials
EGlu = 0;                % Reversal potential of glutamate
EChl = -70;              % Reversal potential of chloride
tau_rise_In = 0.15;      % Tau rise for inhibition
tau_decay_In = 0.9;      % Tau decay for inhibition
tau_rise_Ex = 0.23;      % Tau rise for excitation
tau_decay_Ex = 0.2;      % Tau decay for excitation
tmax = 15;               % Duration of experim

% Set CPSC2
G_max_chl2 = 10;          % Maximal conductance of chloride
G_max_glu2 = 40;          % Maximal conductance of glutamate
Vm2 = -50;           % Array of holding potentials
EGlu2 = 0;                % Reversal potential of glutamate
EChl2 = -70;              % Reversal potential of chloride
tau_rise_In2 = 0.15;      % Tau rise for inhibition
tau_decay_In2 = 0.9;      % Tau decay for inhibition
tau_rise_Ex2 = 0.23;      % Tau rise for excitation
tau_decay_Ex2 = 0.2;      % Tau decay for excitation
tmax2 = 15;               % Duration of experim

% Initialize time 
dt = 0.1;     % time step duration (ms)
t = 0:dt:tmax-dt;  
% Gating variables for inhibition CPSC1
m_plot_I = 1 - exp(-t / tau_rise_In); % Activation
h_plot_I = exp(-t / tau_decay_In); % Inactivation
% Gating variables for inhibition CPSC2
m_plot_I2 = 1 - exp(-t / tau_rise_In2); % Activation
h_plot_I2 = exp(-t / tau_decay_In2); % Inactivation       
% Gating variables for excitation CPSC1
m_plot_E = 1 - exp(-t / tau_rise_Ex); % Activation
h_plot_E = exp(-t / tau_decay_Ex); % Inactivation
% Gating variables for excitation CPSC2
m_plot_E2 = 1 - exp(-t / tau_rise_Ex2); % Activation
h_plot_E2 = exp(-t / tau_decay_Ex2); % Inactivation
% Gating CPSC1
gat_I = m_plot_I .* h_plot_I; % Inhibitory
gat_E = m_plot_E .* h_plot_E; % Excitatory
% Gating CPSC2
gat_I2 = m_plot_I2 .* h_plot_I2; % Inhibitory
gat_E2 = m_plot_E2 .* h_plot_E2; % Excitatory          
% Conductances CPSC1
Gi = G_max_chl .* gat_I; % Inhibitory
Ge = G_max_glu .* gat_E; % Excitatory
% Conductances CPSC2
Gi2 = G_max_chl2 .* gat_I2; % Inhibitory
Ge2 = G_max_glu2 .* gat_E2; % Excitatory

Thank you again for your help.

Sam

[Matt J]

I assume you figured it out (since you Accept-clicked the asnwer)?

[Samuele Bolotta]

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