How to use fmincon to minimize a function with equations and experimental data

Hello,

I have the following problem. I have a set of experimental data and I would like to find the value of two variables which they will minimize the following function:

with the following constraints:

The constants are AE, BE, CE, AW, BW and CW and can be found in the code. The two variables I want to find are the AEW and AWE. So I thought that I can use fmincon with the minimization function in the fun_vl and for the constrains I used the fsolve as it can be solved as a system of non linear equations, if the two variables are given from the fmincon. The experimental data can be found in the following codes:

function y = fun_vl(x)

AEW = x(1);
AWE = x(2);

T_exp = [95.5 89 86.7 85.3 84.1 82.7 82.3 81.5 80.7 79.8 79.7 79.3 78.74 78.41 78.18]; % Celcius
x_exp = [0.019 0.0721 0.0966 0.1238 0.1661 0.2337 0.2608 0.3273 0.3965 0.5079 0.5198 0.5732 0.6763 0.7472 0.8943]; 
y_exp = [0.17 0.3891 0.4375 0.4704 0.5089 0.5445 0.558 0.5826 0.6122 0.6564 0.6599 0.6841 0.7385 0.7815 0.8943];


options.Algorithm = 'levenberg-marquardt';
x0 = ones(1,90);
xopt = fsolve(@root,x0,options);

Topt = xopt(61:75);
yopt = xopt(76:90);

y = sum((5*((Topt - T_exp)/(T_exp))^2) + ((yopt - y_exp)/(y_exp))^2);
function F = root(p)

global AEW
global AWE

AW = 8.07131;
BW = 1730.63;
CW = 233.426;
AE = 8.1122;
BE = 1592.864;
CE = 226.184;
P = 1;

T_exp = [95.5 89 86.7 85.3 84.1 82.7 82.3 81.5 80.7 79.8 79.7 79.3 78.74 78.41 78.18]; % Celcius
x_exp = [0.019 0.0721 0.0966 0.1238 0.1661 0.2337 0.2608 0.3273 0.3965 0.5079 0.5198 0.5732 0.6763 0.7472 0.8943]; 
y_exp = [0.17 0.3891 0.4375 0.4704 0.5089 0.5445 0.558 0.5826 0.6122 0.6564 0.6599 0.6841 0.7385 0.7815 0.8943];

for k = 1:length(T_exp)
    PE(k)     = p(k);
    PW(k)     = p(k+15);
    gammaE(k) = p(k+30);  
    gammaW(k) = p(k+45);  
    T(k)      = p(k+60); 
    y(k)      = p(k+75); 
end
    
for k = 1:length(T_exp)
    
    F(k)   = (1 - y(k))*P - (1 - x_exp(k))*gammaW(k)*PW(k) ;
    F(k+15) = y(k)*P - x_exp(k)*gammaE(k)*PE(k) ;
    F(k+30) = log(gammaE(k)) - AEW*((AWE*(1-x_exp(k)))/(AEW*x_exp(k) + AWE*(1 - x_exp(k))))^2;
    F(k+45) = log(gammaW(k)) - AWE*((AWE*x_exp(k))/(AEW*x_exp(k) + AWE*(1 - x_exp(k))))^2;
    F(k+60) = log10(PE(k)) - (AE - (BE/(CE + T(k))));
    F(k+75) = log10(PW(k)) - (AW - (BW/(CW + T(k))));
end
clc
clear
% Van Laar for ethanol-water

%% Data

T_exp = [95.5 89 86.7 85.3 84.1 82.7 82.3 81.5 80.7 79.8 79.7 79.3 78.74 78.41 78.18]; % Celcius
x_exp = [0.019 0.0721 0.0966 0.1238 0.1661 0.2337 0.2608 0.3273 0.3965 0.5079 0.5198 0.5732 0.6763 0.7472 0.8943]; % 
y_exp = [0.17 0.3891 0.4375 0.4704 0.5089 0.5445 0.558 0.5826 0.6122 0.6564 0.6599 0.6841 0.7385 0.7815 0.8943];

%% Optimization

Nvar = 2;
options = gaoptimset('PopulationSize',50,...
    'FitnessLimit',1e-6,'Generations',50,'TolFun',1e-6,...
    'EliteCount',5,'CrossoverFraction',0.8,'MigrationDirection','both',...
    'MigrationFraction',0.2,'CrossoverFcn',@crossovertwopoint,...
    'FitnessScalingFcn',@fitscalingrank,'PlotFcn',@gaplotbestf);

[P_0,fval]=ga(@fun_vl,Nvar,[],[],[],[],[0 0],[10 10],[]);

options = optimset('Display','iter','PlotFcns',@optimplotfval,'MaxFunEvals',5000);
[P_opt,fval1]=fmincon(@fun_vl,P_0,[],[],[],[],[0 0],[10 10],[],options);

%% Values

disp('-----------------')
disp('Antoine constants')
disp('AEW = ')
disp(P_opt(1))
disp('AWE = ')
disp(P_opt(2))
disp('-----------------')

I get a lot of errors so any help would be great. Basically I want to find the Van Laar coefficients for a mixture of water and ethanol.

[Walter Roberson]

Your use of global is causing problems. You should avoid using global. See http://www.mathworks.com/help/matlab/math/parameterizing-functions.html

Answer 1

[Torsten]

In my opinion, you should reconsider your approach completely.

For suggested values for A_WE and A_EW from the optimizer, you can calculate the 15x1 vector "res" as

res = y_exp - x_exp*gammaE(x_exp)*PE(T_exp) / ((1-x_exp)*gammaW(x_exp)*PW(T_exp) + x_exp*gammaE(x_exp)*PE(T_exp))

and use lsqnonlin to fit A_WE and A_EW so that sum(res.^2) or sum((res/y_exp).^2) becomes minimal.

That's all.

Why do you turn temperature into a free variable in your model ? It's a datum for each of the experiments .