How to use Lagrange multiplier method (maximization) with multiple functions and constraints, for a bioprocess

Hi there, I'm currently trying to build a model to find the optimum conditions for a bioprocess. Extracting glucose and xylose from straw, while keeping the furfural produced below 3 g/l.

The parameters are time (z1), temperature (z2) and concentration (z3) and outputs are furfural, glucose and xylose.

Through genetic programming I have found the following functions:

Furfural: 9.1e-8 exp(-1.0 Time) - 1.33e-4 Concentration Time + 1.01e-6 Concentration Temperature Time + 0.00584

Glucose: (6.78e-26 (6.56e+23 Time - 1.2e+25 Concentration + 1.6e+23 Concentration2 - 9.4e+18 Concentration4 + 1.03e+23 Concentration Temperature))/Concentration

Xylose: 0.0117 Temperature + 0.00229 Time + 0.26 Concentration1/2 - 1.0e-4 Concentration Temperature - 2.0

To find the optimum parameter values I thought Lagrange multiplication method is the way to go.

Constraints: z1, z2, z3 > 0 & furfural production < 3 g/l (this translates to 0.15341 in the formula because the data is in % of aqaeous)

Glucose is twice as "valuable" as xylose, so the maximization function becomes f = 2*glu + xyl

I have found some code online and got to this point, however I'm not sure how to put in the constraints regarding z1, z2 & z3. I don't even know if this code is correct, I'm quite the sucker at coding so please use layman terms :)

My Matlab code looks like this:

syms z1 z2 z3 lambda

%%data

% data in % (not normalized)

glu = (6.781*10^-26*(6.56*10^23*z1 - 1.2*10^25*z3 + 1.6*10^23*z3^2 - 9.4*10^18*z3^4 + 1.03*10^23*z3*z2))/z3 ;% glucose function

xyl = 0.0117*z2 + 0.00229*z1 + 0.26*z3(1/2) - 1.0*10^-4*z3*z2 - 2 ;% xylose function

fur = 9.1*10^-8*exp(-z1) - 1.33*10^-4*z3*z1 + 1.01*10^-6*z1*z2*z3 + 0.00584 + 0.15341 == 0 ; % furfural CONSTRAINT

f = 2*glu + xyl % maximazation function

L = f + lambda * lhs(fur); % Lagrange function

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

system = [dL_dz1; dL_dz2; dL_dz3; dL_dlambda]; % build the system of equations

[z1_val, z2_val, z3_val,lambda_val] = solve(system, [z1 z2 z3 lambda], 'Real', true) % solve the system of equations and display the results

results_numeric = double([z1, z2_val, z3_val,lambda_val]) % show results in a vector of data type double

fur_res = 0.332*z3_val*z2_val^2 - 0.0289*z3_val + 0.332*z2_val*z1_val^2 + 0.332*z1_val*z2_val*z3_val + 0.00825

glu_res = 0.442*z2_val-0.207*z3_val-0.367*exp(-z3_val)*exp(-z1_val)+1.33*z3_val*exp(-z3_val) +0.348

xyl_res = 0.192*z1_val - 0.624*exp(-3*z3_val)+ 0.507*z2_val^(1/2)- 0.11* abs(z3_val)^2*abs(z2_val)* abs(2*z3_val +z1_val) + 0.66

[Max van den Heuvel]

Oh and I am also intersted in having the parameter constraints not z1, z2, z3 > 0 but

10 < z1 < 50

130 < z2 < 170

11 < z3 < 89

Thank you in advance!