How to optimize nonlinear multiple input model with multiple constraints

I'm looking to find optimum parameter values for a bioprocess. Maximise glucose & xylose production, while constraining furfural produced, and make sure glucose & xylose production stay within feasable numbers (function g_glu/g_xyl is in % of aqeaous so it can not be bigger than 1 (100%)).

The model looks something like this.

maximize => f(x1,x2,x3) = 2*f_glu(x1,x2,x3) + f_xyl(x1,x2,x3) (all nonlinear functions)

constraints:

10 < x1 < 50

130 < x2 < 170

11 < x3 < 89

0 < g_glu(x1,x2,x3) < 1

0 < g_xyl(x1,x2,x3) < 1

g_fur(x1,x2,x3) < 3

So far I have tried the live optimisation editor and lagrange but I am not good at writing code and get stuck everytime. I'll dump some of the code that I've writting down below

% LIVE OPTIMISATION EDITOR
function f = objectiveFcn(optimInput)
z1 = optimInput(1);
z2 = optimInput(2);
z3 = optimInput(3);
f = -2*((0.442*z2-0.207*z3-0.367*exp(-z3)*exp(-z1)+1.33*z3*exp(-z3) +0.348) - (0.192*z1 - 0.624*exp(-3*z3)+ 0.507*z2^(1/2)- 0.11* abs(z3)^2*abs(z2)* abs(2*z3 +z1) + 0.66)*(0.916886995261484 - 0.26709662073194) + 0.26709662073194)*1.76 - ((0.192*z1 - 0.624*exp(-3*z3)+ 0.507*z2^(1/2)- 0.11* abs(z3)^2*abs(z2)* abs(2*z3 +z1) + 0.66) * (1.00459544182241 - 0.200345860452685) + 0.200345860452685)*2.958;
end

function [c,ceq] = constraintFcn(optimInput)
% Note, if no inequality constraints, specify c = []
% Note, if no equality constraints, specify ceq = []
z1 = optimInput(1);
z2 = optimInput(2);
z3 = optimInput(3);

c(1) = z1 - 50;
c(2) = 10 - z1;
c(3) = z2 - 170;
c(4) = 130 - z2;
c(5) = z3 - 89;
c(6) = 11 - z3;

ceq = [];
end

% LAGRANGE FORMULA
syms z1 z2 z3 lambda

%conversion rates from normalization to (%)
fur_min = 0.0020923694565838 ;
fur_max = 0.17896582936453 ;
glu_min = 0.26709662073194 ;
glu_max = 0.916886995261484 ;
xyl_min = 0.200345860452685 ;
xyl_max = 1.00459544182241 ;

% functions (normalized)
glu_n = 0.442*z2-0.207*z3-0.367*exp(-z3)*exp(-z1)+1.33*z3*exp(-z3) +0.348 ; % glucose function (normalized)
xyl_n = 0.192*z1 - 0.624*exp(-3*z3)+ 0.507*z2^(1/2)- 0.11* abs(z3)^2*abs(z2)* abs(2*z3 +z1) + 0.66 ; % xylose function (normalized)
fur_n = 0.332*z3*z2^2  - 0.0289*z3  + 0.332*z2*z1^2  + 0.332*z1*z2*z3 + 0.00825 ; % furfural function (normalized)

glu_p = glu_n * (glu_max - glu_min) + glu_min ; % glucose function in (%)
xyl_p = xyl_n * (xyl_max - xyl_min) + xyl_min ; % xylose function in (%)
fur_p = fur_n * (fur_max - fur_min) + fur_min ; % furfural function in (%)

glu = glu_p * 2.958 ; % glucose function in (g)
xyl = xyl_p * 1.584 ; % xylose function in (g)

% constraints
fur = (fur_p * 1.76 / 0.9 * 10)  <= 3 ; % furfural CONSTRAINT in (g/l)
% hb_z1 <= 50
% lb_z1 >= 10
% hb_z2 <= 170
% lb_z2 >= 130          % HOW DO I PUT THESE CONSTRAINTS IN?
% hb_z3 <= 89
% lb_z3 >= 11
% lb_glu >= 0
% hb_glu <= 1
% lb_xyl >= 0
% hb_glu <= 1

% lagrange functions
f = 2*glu + xyl % maximazation function
L = f - lambda * lhs(fur)  % lagrange formula

% calculations
dL_dz1 = diff(L,z1) == 0; % derivative of L with respect to z1 (time)
dL_dz2 = diff(L,z2) == 0; % derivative of L with respect to z2 (temp)
dL_dz3 = diff(L,z3) == 0; % derivative of L with respect to z3 (conc)
dL_dlambda = diff(L,lambda) == 0; % derivative of L with respect to lambda

% outcome
system = [dL_dz1; dL_dz2; dL_dz3; dL_dlambda]; % build the system of equations
[z1_val_n, z2_val_n, z3_val_n,lambda_val] = solve(system, [z1 z2 z3 lambda], 'Real', true)  % solve the system of equations and display the results

Answer 1

[Alan Weiss]

I doubt that you want a symbolic solution. You probably want numbers. So I suggest that you use the Problem-Based Optimization Workflow, using optimization variables instead of symbolic variables.

Alan Weiss

MATLAB mathematical toolbox documentation