Question

Python code for Multivariable Unconstrained *unidirectional search*(topic:- optimization techniques)

Answer 1

Hi,

Please find the code below.

*********************************************************************************

import numpy as np

eps = 0.000001
delta = 0.001
nc = 1 ## number of constraints.
g = np.zeros(nc)

def multi_f(x):
## Himmel Blau function!
sum_ = pow((pow(x[0],2) + x[1] - 11),2) + pow((pow(x[1],2) + x[0] - 7),2)
g[0] = -26.0 + pow((x[0]-5.0), 2) + pow(x[1],2);#constraints.

for i in range(nc):
if(g[i] < 0.0): ## meaning that the constraint is violatd.
sum_ = sum_ + r*g[i]*g[i];

return(sum_)
def grad_multi_f(grad, x_ip):
d1 = np.zeros(M);
d2 = np.zeros(M);
delta_=0.001;
  
for i in range(M):
for j in range(M):
d1[j]=x_ip[j]; d2[j]=x_ip[j];   

d1[i] = d1[i] + delta_;
d2[i] = d2[i] - delta_;

fdiff = multi_f(d1); bdiff = multi_f(d2);   
grad[i] = (fdiff - bdiff)/(2*delta_);
  
return(grad)
## BOUNDING PHASE METHOD,, TO BRACKET THE ALPHASTAR.
def bracketing_(x_0, s_0):
k=0;

## Step__1
## "choose an initial guess and increment(delta) for the Bounding phase method".
del_=0.0000001; w_0=0.5;#initial guess for bounding phase.
  
while(1):
f0 = uni_search(w_0, x_0, s_0, 0.0, 1.0);
fp = uni_search(w_0+del_, x_0, s_0, 0.0, 1.0);
fn = uni_search(w_0-del_, x_0, s_0, 0.0, 1.0);
  
  
if(fn >= f0):
if(f0 >= fp):
del_ = abs(del_); ## step__2
break;
  
else:
a=w_0-del_;
b=w_0+del_;   
  
elif((fn <= f0) & (f0 <= fp)):
del_ = -1*abs(del_);
break;

else:
del_ /= 2.0;

wn = w_0 - del_;## wn = lower limit and w_1 = upper limit.

## Step__3
w_1=w_0 + del_*pow(2,k); ## << exponential purterbation
f1=uni_search(w_1, x_0, s_0, 0.0, 1.0);

## Step__4
while(f1 < f0):
k+=1;
wn=w_0;
fn=f0;
w_0=w_1;
f0=f1;
w_1=w_0 + del_*pow(2,k);
f1=uni_search(w_1, x_0, s_0, 0.0, 1.0);

a=wn;
b=w_1;
## if del_ is + we will bracket the minima in (a, b);
## If del_ was -,
if(b < a):
temp=a;a=b;b=temp;#swapp
return(a, b)
def f_dash(x_0, s_0, xm):
xd = np.zeros(M);

# using central difference.   
for i in range(M):
xd[i] = x_0[i] + (xm+delta)*s_0[i];
fdif = multi_f(xd);

for i in range(M):
xd[i] = x_0[i] + (xm-delta)*s_0[i];
bdif = multi_f(xd);   
  
f_ = (fdif - bdif)/(2*delta);
return(f_)
## This is a function to compute the "z" used in the formula of SECANT Method..
def compute_z(x_0, s_0, x1, x2):
z_ = x2 - ((x2-x1)*f_dash(x_0, s_0, x2))/(f_dash(x_0, s_0, x2)-f_dash(x_0, s_0,x1));
return(z_)

def secant_minima(a, b, x_0, s_0):
## Step__1
x1=a;
x2=b;

##Step__2 --> Compute New point.
z = compute_z(x_0, s_0, x1, x2);
  
##Step__3
while(abs(f_dash(x_0, s_0, z)) > eps):
if(f_dash(x_0, s_0, z) >= 0):
x2=z; # i.e. eliminate(z, x2)   
else:
x1=z; # i.e. elimintae(x1, z)
z = compute_z(x_0, s_0, x1, x2); #Step__2

return(z)
def DFP(x_ip):
eps1 = 0.0001;
eps2 = 0.0001;
  
grad0=np.zeros(M);
grad1=np.zeros(M);
xdiff=np.zeros(M);
Ae=np.zeros(M);
graddiff=np.zeros(M);   
x_1=np.zeros(M);
x_2=np.zeros(M);
  
## step_1
k=0;
x_0=x_ip #This is the initial Guess_.   

## step_2   
s_0=np.zeros(M); ##Search Direction.
grad0 = grad_multi_f(grad0, x_ip);
s_0 = -grad0
  
## step_3 --> unidirectional search along s_0 from x_0. To find alphastar and eventually x(1).
## find alphastar such that f(x + alphastar*s) is minimum.
a,b = bracketing_(x_0, s_0);
alphastar = secant_minima(a, b, x_0, s_0);
  
for i in range(M):
x_1[i] = x_0[i] + alphastar*s_0[i];
  
grad1 = grad_multi_f(grad1, x_1); # Now computing grad(x(1))..

## step_4
## Lets initialize the matrix A,
A = np.eye(M) #Identity Matrix.
dA = np.zeros((M,M))
A1 = np.zeros((M,M))

for j in range(50): # 50 is the total iterations we will perform.   
k+=1;
#print("Iteration Number -- ", k)
for i in range(M):
xdiff[i]=x_1[i]-x_0[i];
graddiff[i]= grad1[i]-grad0[i];
  
# Now, Applying that GIANT FORMULA..
dr=np.inner(xdiff, graddiff);
for i in range(M):
Ae[i]=0;
for j in range(M):
dA[i][j] = (xdiff[i]*xdiff[j])/dr;
Ae[i] += (A[i][j]*graddiff[j]);
  
dr=np.inner(graddiff, Ae)
for i in range(M):
for j in range(M):
A1[i][j] = A[i][j] + dA[i][j] - (Ae[i]*Ae[j])/dr;
A[i][j] =A1[i][j];

## Finding the New Search Direction s_1.. -A_1*grad(x_1)
s_1 =np.zeros(M);
for i in range(M):
s_1[i]=0;
for j in range(M):
s_1[i] += (-A1[i][j]*grad1[j]);

## Unit vector along s_1;
unitv = np.linalg.norm(s_1);
for j in range(M):
s_1[j] /= unitv;
  
## Step__5
# Performing Unidirectional Search along s_1 from x_1; to find x_2;)
a,b=bracketing_(x_1, s_1);
alphastar = secant_minima(a, b, x_1, s_1);
alphastar = (int)(alphastar * 1000.0)/1000.0; #Rounding alphastar to three decimal places.
for i in range(M):
x_0[i] = x_1[i];
x_1[i] = x_1[i] + alphastar*s_1[i];

## Step__6 --->> TERMINATION CONDITION,,
grad1 = grad_multi_f(grad1, x_1);
grad0 = grad_multi_f(grad0, x_0);
  
for j in range(M):
x_2[j] = x_1[j] - x_0[j];
  
if((np.linalg.norm(grad1) <= eps1) or (np.linalg.norm(x_2)/np.linalg.norm(x_0)) <= eps2):
break;

return(x_1)

******************************************************************************

Please contact me if you have any queries.

Please do upvote for my efforts.

Thank you :)