Optimization using Newton Raphson Method

Aim: Determining minimum pressure to break ice using an air cushion ice breaker - Newton Raphson method

Objective:

1) Find minimum pressure for different ice thickness

2) Finding a suitable relaxation factor at h = 0.6 feet

3) Using the optimized relaxation factor to determine minimum pressure and compare it with results obtained with a relaxation factor of 1.

Background:

Air Cushion Vehicle:

In determining the minimum cushion pressure needed to break a given thickness of ice using an air cushion vehicle, Muller (“Ice Breaking with an Air Cushion Vehicle”) in Mathematical Modeling: Classroom Notes in Applied Mathematics, SIAM 1987) derived the equation

`f =p^3(1-beta^2)+p^2(0.4hbeta^2 -(sigmah^2)/r^2)+p((sigma^2h^4)/(3r^4))-((sigmah^2)/(3r^2))^3=0`

`(delf)/(delp) =3p^2(1-beta^2)+2p(0.4hbeta^2 -(sigmah^2)/r^2)+((sigma^2h^4)/(3r^4))`

Numerical solution of equations:

"Numerical methods are often, of a repetitive nature. These consist in repeated execution of the same process where at each step the results of the preceding step is used. This is known as iteration process and is repeated till the result is obtained to a desired degree of accuracy".[1]

Newton Raphson Method: [1]

Let approximate root be `x0` for a function `f(x) = 0`

If `x1 = x0+h` is the exact solution then `f(x1) = 0`

`implies f(x0+h)=0`

Expading using Taylor Series and neglectic ^2 and higher order terms we get:

`f(x0)+hf'(x0) = 0`

`impliesh = - (f(x0))/(f'(x0))`

`implies x1 = x0-(f(x0))/(f'(x0))` Newton-Raphson formula

This is also used with a relaxation factor based on the function to be optimized

`x1 = x0 - rf*(f(x0))/(f'(x0))`

Procedure:

1) Functions 'f' , 'df' are defined using mullers equation (f) and (df/dp). As the pressure is in psi, the function values are of very orders, hence a scaling factor of 1e-26 is used for both functions, therein nullifying the effect when used in numerical method(lines 10-24)

2) Variable t is defined as 7x2 array where column one takes the values of thickness in feet (line 38) and column two carries the corresponding pressure obtained through Newton Raphson method (line 50)

3) Initial approximation of 150 psi is taken and NRM with a relaxation factor of 1 is used to find minimum pressure with a tolerance of 1e-04 for the function f. (lines 34-50)

4) 'time' variable is used to calculate the number of iterations it takes to calculate the min pressures for all values of thickness. (lines 33,49)

5) The results are then plotted and also displayed in a table format (lines 53-64)

6) To find a suitable relaxation factor the thickness is made constant (h = 0.6 feet) And NRM is carried out at different relaxation factors with same tolerance. (lines 71-94)

7) Plots are generated for relaxation factor with respect to number of iterations and the muller function value at minimum pressure. (lines 96-106)

8) Most suitable relaxation factor is the one that can produce the muller function with least deviation from zero (line 89) at the same time arriving at this solution in least number of steps (line 88). To determine the same, the product of both values are compared (lines 90-93) and the minimum value is used to determine the most suitable relaxation factor (lines 108,109)

9) Using the RF obtained above, steps 3-5 have been repeated and number of iterations is stored in time2 (lines 116-144)

10) The difference in number of iterations is calculated and displayed (lines 147,148)

'''
Minimum pressure for air cushion icebreaker
'''

#Library Import
import numpy as np
import matplotlib.pyplot as graph

#function:
def f(p,h,r,beta,sigma):
	a = pow(p,3)*(1-pow(beta,2)) 
	b = ((0.4*h*pow(beta,2))-(sigma*pow(h,2)/pow(r,2)))*pow(p,2)
	c =  pow(sigma,2)*pow(h,4)*p/(3*pow(r,4))
	d = sigma*h*h/3*r*r
	Func_Result = a + b + c - pow(d,3)
	return Func_Result*1e-26

#first derivative of the function
def df(p,h,r,beta,sigma):
	a = 3*p*p*(1-pow(beta,2)) 
	b = 2*p*((0.4*h*pow(beta,2))-(sigma*pow(h,2)/pow(r,2)))
	c = pow(sigma,2)*pow(h,4)/(3*pow(r,4))
	Func_Result = a + b + c
	return Func_Result*1e-26

#Input
beta = 0.5
r = 40 #feet
r = r*12 #inches
sigma = 150 #pound/inch^2(psi)

#pressure minimization for various ice thickness using newton rapson method
time = 0
t = np.zeros((7,2,1))#Initializing table t
for i in range(1,8):
	#thickness array generation
	h = 0.6*i #feet
	t[i-1,0] = h #column 1 of table with thickness in feet
	h = h*12 #converting to inches
	#corresponding minimum pressure array generation
	p = 150 #psi
	j = 1
	tolerance = 1e-04
	while (abs(f(p,h,r,beta,sigma)) > tolerance):
		func = f(p,h,r,beta,sigma)
		dfunc = df(p,h,r,beta,sigma)
		p = p - (func/dfunc)#newton rapsion method
		j = j + 1	
		time = time + 1
	t[i-1,1] = p #colume 2 of table with pressure in psi


#plotting
graph.plot(t[:,0],t[:,1])
graph.plot(t[:,0],t[:,1],'*')
graph.xlabel('Thickness in feet')
graph.ylabel('Pressure in psi')
graph.title('Minimum Pressure vs Thickness')
graph.show()

#printing table
print('|thickness|      pressure      |')
for i in range(0,7):
	print('| ',t[i,0],' | ',t[i,1],' |')



'''
#Relaxation factor optimization at h = 0.6
'''
h = 0.6*12
tolerance = 1e-04
rf = 0.2
relaxationFactor = []
counter = []
funcionVal = []
optim = []
Result = [100,rf]
while (rf < 1.9):
	relaxationFactor.append(rf)
	p = 150
	j = 0
	while (abs(f(p,h,r,beta,sigma)) > tolerance):
		func = f(p,h,r,beta,sigma)
		dfunc = df(p,h,r,beta,sigma)
		p = p - rf*(func/dfunc)
		j = j + 1	
	counter.append(j)
	funcionVal.append(func)
	Op = i*abs(func)
	optim.append(Op)
	if (Result[0] > Op):
		Result = [Op,rf]
	rf = rf + 0.05

graph.plot(relaxationFactor,counter)
graph.plot(relaxationFactor,counter,'o')
graph.xlabel('Relaxation Factor')
graph.ylabel('Number of Iterations')
graph.show()

graph.plot(relaxationFactor,funcionVal,'x')
graph.plot([0,2],[0,0])
graph.xlabel('Relaxation Factor')
graph.ylabel('Tolerance achieved')
graph.show()

print('Optimized Relaxation Factor:')
print(Result[1])



'''
pressure minimization for various ice thickness using newton rapson method with optimized relaxation factor
'''
ORF = Result[1]
time2 = 0
t2 = np.zeros((7,2,1))#Initializing table t
for i in range(1,8):
	#thickness array generation
	h = 0.6*i #feet
	t2[i-1,0] = h #column 1 of table with thickness in feet
	h = h*12 #converting to inches
	#corresponding minimum pressure array generation
	p = 150 #psi
	j = 1
	tolerance = 1e-04
	while (abs(f(p,h,r,beta,sigma)) > tolerance):
		func = f(p,h,r,beta,sigma)
		dfunc = df(p,h,r,beta,sigma)
		p = p - ORF*(func/dfunc)#newton rapsion method
		j = j + 1	
		time2 = time2 + 1
	t2[i-1,1] = p #colume 2 of table with pressure in psi


#plotting
graph.plot(t2[:,0],t2[:,1],'b-')
graph.plot(t[:,0],t[:,1],'r--')
graph.xlabel('Thickness in feet')
graph.ylabel('Pressure in psi')
graph.legend(['Optimized relaxation','Relaxation = 1'])
graph.title('Minimum Pressure vs Thickness with optimized relaxation factor')
graph.show()


print('Time saved using relaxation factor in terms of number of iterations: ')
print(time - time2)

Results & Discussion:

1)Plot 1:

2)Plot 2:

3)Plot 3:

4)Plot 4:

5)Command Window:

The plots and results show that by using a suitable relaxation factor we can obtain results with as much accuracy, in a shorter processing time. We can also observe that the RF is not directly related to the number of iteration. In fact the relation is highly dependent on the function being optimized.

References:

1) Higher Engineering Mathematics by B.S.Grewal - 44th Edition