What is wrong with this optimization code using fmincon?

1 view (last 30 days)
Diego Torres-Siclait on 23 Feb 2021
Answered: Walter Roberson on 23 Feb 2021
Baisically using fmincon to solve for an unknown 2x1 matrix called A, which is optimized value for given objective function. Code error appears in 1st line below but I want A to be a 2x1. For context, x, y, and P are all vectors of the same length.
% Model section
GeRT_model = @(A) (x.*(1.-x)).*[(1.-x) x]*A;
% Parameters for experimental section below. Lots of distracting numbers.
A1 = 7.89750; B1 = 1474.08; C1 = 229.13;
A2 = 8.01195; B2 = 1698.785; C2 = 231.04;
T = 333.15;
P1sat = 10^( A1 - B1 / ((T-273.15) + C1) ) ;
P2sat = 10^( A2 - B2 / ((T-273.15) + C2) );
gamma1 = y.*P ./ (x.*P1sat);
gamma2 = (1.-y).*P ./ ( (1-x).*P2sat );
% experimental section
GeRT_exp = x.*log(gamma1) + (1.-x).*log(gamma2);
%% Optimization: Objective Function that is defined in terms of model and exp. sections. N is length of P,x,and y vectors. The rest is to set up Optimization
OBJFUN = @(A) 1/N *sum ( ( ( GeRT_model(A) .- GeRT_exp ) ./ GeRT_exp ).^2 );
% ineq constrnts
A = [];
b = [];
% eq constrnts
%[x2^2*(1-2*x1) 2*x1*x2^2; 2*x1^2*x2 x1^2*(1-2*x1)]*...
%[A12 A21]'=[log(gamma1) log(gamma2)]'
Aeq = [(1.-x).^2.*(1.-2.*x) 2.*x.*(1.-x).^2; 2.*x.^2.*(1.-x) x.^2.*(1.-2.*x)];
beq = [log(gamma1) log(gamma2)]';
lb = [];
ub = [];
A0 = [1 1]; % an initial guess for A
[A,fval] = fmincon(OBJFUN,A0,A,b,Aeq,beq,lb,ub)
Any help is much appreciated!!

Walter Roberson on 23 Feb 2021
OBJFUN = @(A) 1/N *sum ( ( ( GeRT_model(A) .- GeRT_exp ) ./ GeRT_exp ).^2 );
^^
.- is not a valid operation.
A0 = [1 1]; % an initial guess for A
That is 1 x 2. That is providing an initial value for the A that is being passed into OBJFUN
OBJFUN = @(A) 1/N *sum ( ( ( GeRT_model(A) - GeRT_exp ) ./ GeRT_exp ).^2 );
OBJFUN passes it into GeRT_model
GeRT_model = @(A) (x.*(1.-x)).*[(1.-x) x]*A;
GeRT_model takes it and does matrix multiplication of something with A. With A being 1 x 2, that means the first value must be N x 1 for some N, giving an N x 2 result. But we can see several cases:
• if x is empty, then [(1.-x) x] is empty and the * operation is invalid
• if x is a column vector then [(1.-x) x] is N x 2 and the * operation against 1 x 2 A is invalid
• if x is a row vector then [(1.-x) x] is 1 x (2*N) and the * operation against 1 x 2 A is invalid
so we must figure that you mean
GeRT_model = @(A) (x.*(1.-x)).*[(1.-x) x]*A.';
and that you expect x to be a column vector.
Aeq = [(1.-x).^2.*(1.-2.*x) 2.*x.*(1.-x).^2; 2.*x.^2.*(1.-x) x.^2.*(1.-2.*x)];
beq = [log(gamma1) log(gamma2)]';
Those constraints are unlikely to be satisfied. A will be a vector of two values, and Aeq*A' = beq must be true to within tolerance for all 2*N entries. Your x, y, P would have to have magic relationships with essentially no noise.