Reputation: 11
I'm trying to numerically solve some mechanism design question using the following code. However, the output is not optimal(just stay at my initial guess). How should I make it work?
I tried different initial guesses ( I can calculate myself so the initial guess is valid, each x0 may give different solutions).
I also switched methods, but some results (like Nelder-Mead) generate really invalid solutions (violating my constraints).
import numpy as np
from scipy.optimize import minimize
p0=0.3
r1=0.4
r2=0.8
p1s=r1/(r1+1)
p2s=r2/(r2+1)
q1=(r1+r1*r2)/(r2+r1*r2)
q2=r1/r2
q=(q1+q2)/2
def v_1(k,p0):
if p0>p1s and k>p0:
return r1-p0*r1/k
elif p0<p1s and k>p1s:
return p0+p0*r1-p0*r1/k
return 0
def v_2(k,p0):
if p0>p2s and k>p0:
return r2-p0*r2/k
elif p0<p2s and k>p2s:
return p0+p0*r2-p0*r2/k
return 0
def V(p):
if p<p1s:
return p
if p>p2s:
return q*r2*(1-p)
A=r1-r1/(r2-r1)*(p+p*r2-r1)
return q*p+(1-q)*A
def objective(x):
k1,k2,p1,p2 = x
return -q*(p2+V(k2))+(1-q)*(p1+V(k1)) # Convert to minimization
def IR1(x):
k1,k2,p1,p2 = x
return v_1(k1,p0)-p1
def IR2(x):
k1,k2,p1,p2 = x
return v_2(k2,p0)-p2
def IC1(x):
k1,k2,p1,p2 = x
return v_1(k1,p0)-p1-v_1(k2,p0)-p2
def IC2(x):
k1,k2,p1,p2 = x
return v_2(k2,p0)-p2-v_2(k1,p0)-p1
constraints = [
{'type': 'ineq', 'fun': lambda x: IR1(x)},
{'type': 'ineq', 'fun': lambda x: IR2(x)},
{'type': 'ineq', 'fun': lambda x: IC1(x)},
{'type': 'ineq', 'fun': lambda x: IC2(x)}
]
# Bounds for x1, x2
bounds = [(p0, 1), (p0, 1), (0,p0), (0,p0)] #k1,k2>=p0,p0>=p1,p2>=0
# Initial guess
x0 = [p0, 1, 0,p0 ]
# Solve the problem
result = minimize(objective, x0, bounds=bounds, constraints=constraints,method='SLSQP')
# Print results
if result.success:
print("Optimal solution found:")
print(f"k1 = {result.x[0]:.4f}, k2 = {result.x[1]:.4f}")
print(f"p1 = {result.x[2]:.4f}, p2 = {result.x[3]:.4f}")
print(f"Maximum function value: {-result.fun:.4f}")
else:
print("Optimization failed:", result.message)**strong text**
Upvotes: 0
Views: 22