import math
from pfdLib import PFD

c_I =	10000
c_H =	10000000
c_R =	50000
lmbda =	0.01752
f_D =	0.2
beta =	0.1
k=1
n=2


def C(tau,lmbda,f_D,c_I,c_R,c_H,k,n,beta):
    return c_I/tau+n*c_R*(lmbda-lmbda**2*tau/2)+c_H*PFD(tau,lmbda,k,n,beta)*f_D

from scipy.optimize import minimize    
tau0 = 1.1 # Initial guess

result = minimize(C, tau0, args = (lmbda,f_D,c_I,c_R,c_H,1,1,0))
print("\nTest for k=1, n=1")
print("tau* =",result.x[0])
print("C(tau*) =",result.fun)


 
result = minimize(C, tau0, args = (lmbda,f_D,c_I,c_R,c_H,k,n,beta))
print("\nTest for k=1, n=2")
print("tau* =",result.x[0])
print("C(tau*) =",result.fun+15000)