Matlab Program for Quadratic Interpolation

“Quadratic Interpolation” program in  MATLAB for one-dimensional minimization. Program is a generalized one and takes objective function as user input in terms of standard MATLAB function, e.g. for minimization of f(x)= 3×2-8x+2, user will enter ‘3*x^2-8*x+2’. For theoretical understanding of Quadratic Interpolation you can follow: Linear and quadratic interpolation

The inputs to the program are as follows:

S. No. Input Remarks
1 objective function in string form e.g. ‘3*x^2-8*x+2’
2 Lower Limit of x
3 Upper Limit of x
4 Initial step size
5 Max Iteration count for inner loop Reduction in step length for learning of a b and c
6 Max Iteration count for outer loop Maximum number of times a b and c can be searched
7 Min value of step length t0
8 Epsilon value for stopping criterion 1 (h-f)/f
9 Epsilon value for stopping criterion 2 takes difference of left and right hand side derivative

Code: (code is free to use with acknowledgement)

function [ ] = QuadInterpolation( )
%Quadratic Interpolation method 
% Quadratic Interpolation method for one dimensional minimization problem
% user Inputs are:
% 1) objective function in string form (e.g. 3*x^2-8*x+2)
% 2) function domain(lower and upper limit of x)
% 3) initial step size to
% 4) Maximum Iteration count before reinitializing values of a ,b and c
% 5) Maximum number of times a,b, and c can be reinitialized
% 6) Epsilon for relative value stopping criterion
% 7) Epsilon for derivative based stopping criterion
% Written By: Himanshu 
clc;

% start segment user inputs and equation parsing-----------------------------------------
enteryString='Generalized code for Quadratic Interpolation in one-dimension. \n Enter objective Function (example 3*x^2-8*x+2):\n';
obfunStr=input(enteryString);
obFunEval=sym(obfunStr);
Llimit=input('Enter Lower Limit of x: ');
Ulimit=input('Enter Upper Limit of x: ');
t0=input('Enter Initial Step size: ');
N=input('Enter Maximum Iteration count before reinitializing values of a ,b and c: ');
M=input('Enter Maximum number of times a,b, and c can be reinitialized: ');
minStep=input('Enter minimum value of step Length: ');
epsRel=input('Enter Epsilon for "relative value based" stopping criterion: '); % (h-f)/f
epsDer=input('Enter Epsilon for "derivative based" stopping criterion: ');%(fdellamL-fdellamR)/2) takes average of left and right hand side derivative as stopping criterion

% end segment user inputs and equation parsing-------------------------------------

bestVals=nan(M,2); % stores best value and corresponding x

%plot objective function for min 100 sampling points------------------------
samPoints=max(ceil(Ulimit-Llimit),100)+1;
yVals=nan(samPoints,1);
xVals=nan(samPoints,1);
for i=1: samPoints 
 xValss=(Llimit+((Ulimit-Llimit)*(i-1)/(samPoints-1)));
 if(xValss~=0)
 yVals(i)=subs(obFunEval,'x',(Llimit+((Ulimit-Llimit)*(i-1)/(samPoints-1))));
 xVals(i)=(Llimit+((Ulimit-Llimit)*(i-1)/(samPoints-1)));
 end
 
end
plot(xVals,yVals,'r','LineWidth',2)
xlabel('x');
ylabel('f(x)');
title(strcat('Optimization function : ',obfunStr));
hold on;
grid;
%--------------------------------------------------------------------------

for i=1:M
 %generate the value of a(starting point) randomly between Llimit and Ulimit
 disp(strcat('Iteration Count (Outer Loop): ',num2str(i)));
 a= Llimit + (Ulimit-Llimit).*rand(1,1);
 % code segment for learning values of a b and c------------------------
 b=nan;
 c=nan;
 stepReversed=false; % set value of stepReversed to false. it is set to true when t0 takes negative value
 pointsFound=true; % this will be set to false when even with reversing t0 points are not found
 initialPoint=a;
 stepLength=t0;
 newPoint=a+stepLength;
 if(subs(obFunEval,'x',initialPoint)<=subs(obFunEval,'x',newPoint)) %faliure
 while (subs(obFunEval,'x',initialPoint)<=subs(obFunEval,'x',newPoint))%reduce step length
 stepLength=stepLength/2;
 if ((stepLength<minStep) && stepReversed==false)% reduce till step length is less than minStep otherwise reverse step
 stepReversed=true;
 stepLength=-t0;
 end
 if ((abs(stepLength)<minStep) && stepReversed==true)% if by reversing step not moving further 
 pointsFound=false;
 break; % stop while loop
 end
 newPoint=initialPoint+stepLength;
 end
 if (pointsFound==false)
 continue; % skip current iteration and find new values for a b and c
 end
 if (stepReversed==true) % reorder a b and c in order a<b<c if t0 is set to negative
 c=initialPoint;
 b=newPoint;
 a=b+stepLength;
 else
 b=newPoint;
 c=b+stepLength;
 a=initialPoint;
 end
 else % success
 temp=initialPoint;
 while (subs(obFunEval,'x',temp)>subs(obFunEval,'x',newPoint))%increase step length
 stepLength=stepLength*2;
 temp=newPoint;
 newPoint=initialPoint+stepLength;
 end
 a=initialPoint;
 b=temp;
 c=newPoint;
 end
 
 % if values of c is greater than upper limit make it equal to upper limit
 if (c>Ulimit) 
 c=Ulimit;
 end
 % if values of b is greater than upper limit search for new a b c
 if (b>Ulimit)
 continue;
 end
 
 % if values of c is less than lower limit make it equal to lower limit
 if (c<Llimit)
 c=Llimit;
 end
 
 % if values of b is less than lower limit search for new a b c
 if (b<Llimit)
 continue;
 end
 
 % display learned values of a b and c
 disp(strcat('Learned Values of a b c for iteration count: ',num2str(i),'a=',num2str(a),' b=',num2str(b),' c=',num2str(c)));
 
 lambda=nan;
 for j=1:N % inner loop one a b and c are determined shring the interval
 funMat=vpa([vpa(subs(obFunEval,'x',a));vpa(subs(obFunEval,'x',b));vpa(subs(obFunEval,'x',c))]);
 % calculate A B and C
 QuadFitParams=vpa(([1 a a^2;1 b b^2;1 c c^2])\funMat);
 lambda=-QuadFitParams(2)/(2*QuadFitParams(3));
 if (QuadFitParams(3)<0) % if c<0 find new values of a b and c
 break; 
 end
 
 %test for covergence
 %calculate h(lambda)
 h=QuadFitParams(1)+QuadFitParams(2)*lambda+QuadFitParams(3)*lambda^2;
 %calculate f(lambda)
 f=vpa(subs(obFunEval,'x',lambda));
 fdellamL=vpa(subs(obFunEval,'x',(lambda+0.01))); 
 fdellamR=vpa(subs(obFunEval,'x',(lambda-0.01)));
 %test convergence
 if(abs(((h-f)/f))<epsRel && abs((fdellamL-fdellamR)/2)<epsDer) 
 disp('stop criterion hit'); % step out of loop
 break; 
 else % shrink interval size by reassigning a b and c
 if (lambda<b && f>funMat(2))
 a=lambda;
 end
 
 if (lambda<b && f<funMat(2))
 c=b;
 b=lambda;
 end
 
 if (lambda>b && f>funMat(2))
 c=lambda;
 end
 
 if (lambda>b && f<funMat(2))
 a=b;
 b=lambda;
 end
 end
 end
 bestVals(i,1)=vpa(subs(obFunEval,'x',lambda));
 bestVals(i,2)=lambda;
 
 
end
% display optimization results
disp('Optimization Results: ');
FinalResults=nan;
for j=1:size(bestVals,1)
if(~isnan(bestVals(j,1)) && ~isnan(bestVals(j,2)) && bestVals(j,2)>=Llimit && bestVals(j,2)<=Ulimit)
FinalResults(j,1)=bestVals(j,1);
FinalResults(j,2)=bestVals(j,2);
disp(strcat('x* = ',num2str(FinalResults(j,2)),' f(x*)=',num2str(FinalResults(j,1))));
% plot x* and f(x*) for current set of a b and c
plot(FinalResults(j,2),FinalResults(j,1),'o','MarkerEdgeColor','k','MarkerFaceColor','g',...
 'MarkerSize',10);
end
end

[bestVal, indexbestval]=min(FinalResults(:,1));
% plot best x* and f(x*) from different value of a b and c
plot(FinalResults(indexbestval,2),bestVal,'s','MarkerEdgeColor','k','MarkerFaceColor','m',...
 'MarkerSize',10);
 end


Lets Evaluate Program on a test function: (x2-5x+6)/(x2+1)

XL: 0

XU: 2*pi

t0:0.5

N:10

M:5

tmin:0. 05

Eps1:0.001

Eps2:0.001

x* =2.4489 f(x*)=-0.035356

 

optimresults

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s