# How to limit the output value of a function to positive values only?

11 views (last 30 days)
Dear MATLAB community,
I am workig on a non-linear regression problem. I already coded my script and the optimization works as well. However, I am facing a problem here in terms of the parameters that GlobalSearch solver sets as optimal solutions. Reason being, they are not nicely set to work well with the model.
I am using the GlobalSearch-solver for determining optimal fittings parameters so that the difference between my model prediction and the original value obtained from experiments is as close to zero as possible.
My model is a cross-model that aims at predicting a viscosity in dependence of a shear rate gamma_dot. It is composed out of two functions: eta0_pr and etainf_pr. BOTH of these functions are only allowed to return positive values. Whereas the output of eta0_pr must always be smaller than eta_inf in any case. This is the shape of the model:
((etainf_pr+((eta0_pr)-(etainf_pr))./(1+x(48).*gamma_dot.^x(49))))
The problem that I am facing is that through the optimization by the solver, my fitting parameters are set in a way that
2. work against mathematical logic. Hence, my 'Initial value' and 'objective value' might be negative. Even though, my objective function is the sum of a vector output which should in any case be positive because of the ^2.
How is this possible?
Is there a way to limit the output of a function in a way that is similar to the MAX-function of Excel? Example:
MAX(f(x),g(x)) and then something like MAX(g(x),c)
where f(x) and g(x) are functions and c is a scalar.
I probably need something like this to say that both 'eta0_pr' and 'etainf_pr' must be bigger or equal to zero PLUS 'eta0_pr' must be bigger than 'etainf_pr'
This is my code:
%% creating an objective function
eta0_pr = @(x,cststp,cst,coil,cegg,rpm,ctsalt,csuc,pH)((((x(1)+x(2).*(cststp-x(3)).^(x(4)).*(cst./100)))+...
(x(5)+x(6).*(coil-x(7)))+(x(8).*(pH-x(9)).^(-1))+(x(10).*csuc)+(x(11).*ctsalt)+...
(x(12).*(cegg-x(13))-(x(14).*(rpm-x(15)))+(x(16).*exp((x(17)+x(18).*(coil-x(7)))+...
((x(19)+x(20).*(cststp-x(3).^x(21))).*(cst./100)))+...
(cegg.*x(22).*(pH-x(23)).^(-1))+(ctsalt.*x(24).*(pH-x(25)).^(-1))+...
(cegg.*ctsalt.*x(26).*(pH-x(27).^(-1)))+(coil.*x(28).*(pH-x(29)).^(-1))))));
etainf_pr = @(x,cststp,cst,coil,rpm,ctsalt,csuc,pH,eta_h2o)((x(30)+x(31).*(cststp-x(3)).^x(32)).*(cst./100))+...
(x(33)+x(34).*(coil-x(35)))-(x(36).*(rpm-x(37)))+(x(38).*ctsalt)+(x(39).*csuc)+...
(x(40).*exp((x(41)+x(42).*(coil-x(35)))+((x(43)+x(44).*(cststp-x(3)).^x(45)).*(cst/100))))+...
(x(46).*(pH-x(47)).^(-1))+(eta_h2o);
eta_pr =@(etainf_pr,eta0_pr,x,gamma_dot)((etainf_pr+((eta0_pr)-(etainf_pr))./(1+x(48).*gamma_dot.^x(49))));
variance =@(eta_pr,eta_m) ((((eta_pr-eta_m))).^2);
objective=@(variance) sum(variance);
% defining initial fitting parameters
x0 = initial_parameters(2, :);
% display initial objective value
disp(['initial objective: ' num2str(objective(x0))])
%% defining problem structure for GlobalSearch solving method
A = [];
b = [];
Aeq = [];
beq = [];
lb = initial_parameters(3, :);
ub = initial_parameters(4, :);
nonlcon =[];
opts = optimoptions(@fmincon,'Display','iter','Algorithm','interior-point'); % choosing solver algorithm
%% fmincon test
x = fmincon(objective,x0,A,b,Aeq,beq,lb,ub,nonlcon,opts);
disp(x)
disp(['final objective: ' num2str(objective(x))])
%% creating problem structure for GlobalSearch solver to run with
problem = createOptimProblem('fmincon','x0',x0,'objective',objective,'Aineq',A,'bineq',b,'Aeq',Aeq,'beq',beq,'ub',ub,'lb',lb,'nonlcon',nonlcon,'options',opts);
%% settings for Global Search solver
gs = GlobalSearch('NumTrialPoints',10000,'NumStageOnePoints',2000,'FunctionTolerance',1.0000e-8,'XTolerance',1.0e-07);
gs
%% run Global Seearch solver
[x,f,exitflag,output,solutions] = run(gs,problem);
%% display resutls
disp(x)
disp(['final objective: ' num2str(objective(x))])
I also uploaded the data and the full code, so if you wanted you could run the script yourself. May be you easily spot the mistake I am making here.
As you can see from the code, I set upper and lower bounds. So, that fmincon has a certain range for all fitting parameters to reiterate. From my understanding, when using GlobalSearch on an objective function it will look for the lowest minimum possible of the specified function by, in this case, testing out different values for 'x' (fitting parameters).
What confuses me is the fact that the output of the objective-function turns negative. Even though, when I compute the same in Excel it does not. A sum of something positive should also always be positive.
Has it to do something with complex numbers in my calculations? that would be the only possibility why a squared value is still negative.
Thanks in advance for any input of yours!
p.s. I am sorry for the clusterf*ck of text in eta0_pr and etainf_pr. It's definitely not reader-friendly

Alan Weiss on 23 Apr 2020
It is difficult to understand what you mean when you say "they are not nicely set to work well with the model." I also am not sure why you use GlobalSearch for solving a least squares problem.
For an example of using MultiStart to find better local solutions, see MultiStart using lsqcurvefit or lsqnonlin.
Even if you really need to use fmincon as your local solver, I suggest that you try MultiStart as the global solver.
Alan Weiss
MATLAB mathematical toolbox documentation
Okay, I will try using MultiStart