python - 使用 scipy.optimize 的 Python 约束非线性优化无法找到最佳解决方案

Question

我有一个问题，我试图找到一个解决方案，因为我有 5 个单变量多项式，这些多项式在我关心的范围内有一个峰值。我的目标是找到每个多项式的变量值（在某些最小/最大值和所有变量约束的总和下），使这些曲线的值乘以每条曲线的常数。

我已经使用 scipy.optimize 包和 numpy 设置了一些代码。它似乎能够达到解决方案，但它所达到的解决方案似乎并不接近最优。例如，普通情况是输入 488 MW。这个特定的输入值有一个解决方案，其中每个 x0-x4 变量都处于其函数的峰值，如下所示：

x0=90
x1=100
x2=93
x3=93
x4=112

The result it provides we with is:

x0=80
x1=97
x2=105
x3=80
x4=126

This does satisfy my constraint, but it does not appear to minimize the objective function. 

import numpy as np
import matplotlib.pyplot as plt

from scipy.optimize import minimize



U1_Max=103
U1_Min=80
U2_Max=102
U2_Min=80
U3_Max=105
U3_Min=80
U4_Max=100
U4_Min=80
U5_Max=126
U5_Min=90


# Our whole goal here is to maximze the sum of several efficiency efficiency curves times the
# output MW of each unit. "The most efficiency where it matters the most"
# Assuming all units are available for assignment his would look something like:

#Where have the following efficiency curves:
#U1: Efficiency=-0.0231*(MW)^2+4.189*MW-102.39
#U2: Efficiency= -0.01*(MW)^2+1.978*MW-8.7451
#U3: Efficiency= -0.025*MW^2+4.5017*MW-115.37
#U4: Efficiency= -0.01*(MW)^2+1.978*MW-8.7451
#U5: Efficiency= -0.0005*(MW)^2+0.1395*(MW)^2-13.327*MW+503.41

#So I think we want to
#Maximize U1(x0)*U1_MAX+U2(x1)*U2_MAX+U3(x2)*U3_MAX+U4(x3)*U4_MAX+U5(x4)*U5_MAX
#I think this can also be stated as:
#Minimize (U1(x0)-100)*U1_MAX+(U2(x1)-100)*U2_MAX)+(U3(x2)-100)*U3_MAX)+(U4(x3)-100)*U4_MAX)+(U5(x4)-100)*U5_MAX)
#Which means 'minimize the sum of the products of the difference between 100% efficient and actual and the unit nameplates'

#By Choosing {x1, x2, x3, x4, x5}

#Such that x1+x2+x3+x4+x5=MW_Total
#Such that U1_Min<x1<U1Max
#Such that U2_Min<x2<U2Max
#Such that U3_Min<x3<U3Max
#Such that U4_Min<x4<U4Max
#Such that U5_Min<x5<U5Max

##so let's type that out....
#stack overflow says the optimizer does best if the function being optimized is around 1-5ish so we will get it there-ish. 
def objective(x):
  return (
      (
      ((100-0.0231*x[0]**2+4.189*x[0]-102.39))*U1_Max
      +((100-0.01*x[1]**2+1.978*x[1]-8.7451))*U2_Max
      +((100-0.025*x[2]**2+4.5017*x[2]-115.37))*U3_Max
      +((100-0.01*x[3]**2+1.978*x[3]-8.7451))*U4_Max
      +((100-0.0005*x[4]**3+0.1395*x[4]**2-13.327*x[4]+503.41))*U5_Max
      )

      )

x=np.zeros(5)

print(
      (
      ((100-0.0231*x[0]**2+4.189*x[0]-102.39))*U1_Max
      +((100-0.01*x[1]**2+1.978*x[1]-8.7451))*U2_Max
      +((100-0.025*x[2]**2+4.5017*x[2]-115.37))*U3_Max
      +((100-0.01*x[3]**2+1.978*x[3]-8.7451))*U4_Max
      +((100-0.0005*x[4]**3+0.1395*x[4]**2-13.327*x[4]+503.41))*U5_Max
      )

      )

#Now, let's formally define our constraints
#Note that this must be of a form that satisfies 'constraint equal to zero'
#First, the sum of all MW commands should be qual to the total MW commanded
def constraint1(x):
    return -x[0]-x[1]-x[2]-x[3]-x[4]+MW_Total

#Since this is a numeric process let's give it some starting 'guess' conditions.
n=5
x0=np.zeros(n)
x0[0]=90
x0[1]=100
x0[2]=93
x0[3]=93
x0[4]=112

# show initial starting uess
print('Start by guessing: ')
print(x0)
print('Which gives a scaled algorithim value of: ')
print(
      (
      ((100-0.0231*x0[0]**2+4.189*x0[0]-102.39))*U1_Max
      +((100-0.01*x0[1]**2+1.978*x0[1]-8.7451))*U2_Max
      +((100-0.025*x0[2]**2+4.5017*x0[2]-115.37))*U3_Max
      +((100-0.01*x0[3]**2+1.978*x0[3]-8.7451))*U4_Max
      +((100-0.0005*x0[4]**3+0.1395*x0[4]**2-13.327*x0[4]+503.41))*U5_Max
      )

      )
print('Which gives actual MW total of: ')
print(x0[0]+x0[1]+x0[2]+x0[3]+x0[4])


#Next, Let's give it some bounds to operate in
U1_Bnds=(U1_Min, U1_Max)
U2_Bnds=(U2_Min, U2_Max)
U3_Bnds=(U3_Min, U3_Max)
U4_Bnds=(U4_Min, U4_Max)
U5_Bnds=(U5_Min, U5_Max)
Bnds=(U1_Bnds, U2_Bnds, U3_Bnds, U4_Bnds, U5_Bnds)

con1 = {'type': 'eq', 'fun': constraint1}
print('MW Generated is: ')
for i in range (410,536):
  MW_Total=i
  solution = minimize(objective,x0,method='SLSQP',bounds=Bnds,constraints=con1,options={'maxiter': 10000000, 'eps':1.4901161193847656e-10})
  x = solution.x
  print(solution.x[0],solution.x[1],solution.x[2],solution.x[3],solution.x[4])

我希望对于我的 488 兆瓦的琐碎案例，它会给我最佳答案。我究竟做错了什么？

标签： pythonnumpyscipymathematical-optimization

Answer 1

通过查看您的目标和约束定义，您似乎处于具有线性约束的二次目标函数的情况。

这方面的理论是众所周知的，并提供了收敛保证，您可以参考维基百科页面。

我不太了解 scipy SLSQP 界面，但看起来你使用的信息比你能做的要少。尝试以具有线性约束的二次目标函数的形式来解决您的问题。还将您的约束投射到一个scipy.optimize.LinearConstraint对象中。

并且请在您的代码中使用诸如print(objective(x))和之类的函数调用print(solution.x)，这将增强可读性。