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** Continue reading...
