% clear session, close plots, clear screen
clear all; close all; clc
% data for regression
xm = [18.3447,79.86538,85.09788,10.5211,44.4556, ...
               69.567,8.960,86.197,66.857,16.875, ...
               52.2697,93.917,24.35,5.118,25.126, ...
               34.037,61.4445,42.704,39.531,29.988];
ym = [5.072,7.1588,7.263,4.255,6.282, ...
               6.9118,4.044,7.2595,6.898,4.8744, ...
               6.5179,7.3434,5.4316,3.38,5.464, ...
               5.90,6.80,6.193,6.070,5.737];
% initial parameter guess
p0 = [1,0.5,1];           
% define prediction function
yp = @(p) p(1) + p(2)./xm + p(3).*log(xm);           
% define objective function (scaled sum of squared errors)
objective = @(p) sum(((yp(p)-ym)./ym).^2);
disp(['Initial Objective: ' num2str(objective(p0))])

% optimize with fmincon
%[X,FVAL,EXITFLAG,OUTPUT,LAMBDA,GRAD,HESSIAN] 
% = fmincon(FUN,X0,A,B,Aeq,Beq,LB,UB,NONLCON,OPTIONS)
A = [];
b = [];
Aeq = [];
beq = [];
% bounds
lb = []; % ones(3)*0.2;
ub = []; % ones(3)*1.5;
popt = fmincon(objective,p0,A,b,Aeq,beq,lb,ub);

% print results
disp(['Final Objective: ' num2str(objective(popt))])
disp(['Optimal parameters: ' num2str(popt)])

% plot results
plot(xm,ym,'ro')
hold on
plot(xm,yp(p0),'bx')
plot(xm,yp(popt),'gs')
legend('measured','initial predicted','optimal predicted')
ylabel('y')
xlabel('x')