首页 > 编程问答 >无法在 Gekko 中求解 MINLP（警告不再有可能的试验点且没有整数解）

无法在 Gekko 中求解 MINLP（警告不再有可能的试验点且没有整数解）

时间：2024-08-07 03:49:31浏览次数：14  
标签：python gekko

我对 python 中的 Gekko 包非常陌生。我的目标是最大化 'Q_factor' 我已有的训练有素的 TensorFlow 模型 (.keras)。

我在稳态模式 ( m.options.IMODE=2 ) 下使用 Gekko 进行参数估计，并使用 ( m.options.SOLVER=1 ) 进行以下一些约束：

我的 MINLP 问题

注意 N 必须是 2 到 8 之间的整数。

1. 首先，我可以解决上面的问题*但结果 -m.options.OBJFCNVAL 看起来不正确（解决方案太大）。解决方案 pred_Q_factor 应该在~700-1000左右。你可以看到下面的代码和解决方案： 

## 1st CODE ##

import numpy as np
import pandas as pd
from gekko import GEKKO
from gekko.ML import Gekko_NN_TF, CustomMinMaxGekkoScaler
from tensorflow.keras.models import load_model
import warnings
warnings.filterwarnings("ignore")

# Q-factor prediction from trained TF model
model_path = 'model_path'
model = load_model(model_path)

# CSV file reading
data_path = 'datapath'
data = pd.read_csv(data_path)

# Define output feature of the model
target = ['Q_factor']
# Basically, the input features are [Freq, Do_1, Do_2, d, p1, p2, N, w1, w2, Di_1, Di_2, t1, t2]
features = [col for col in data.columns if col not in target]

# Find max/min of scaler for Gekko_NN_TF parameters
s = CustomMinMaxGekkoScaler(data, features, target)
mma = s.minMaxValues()

# Initialize GEKKO
m = GEKKO(remote=False)

# Design Requirements
freq_const = 13.56
Do_const = 200
sys_thickness = 200
copper_thickness = 0.5

# Define variables with bounds
Freq = m.FV(value=freq_const, lb=1, ub=30,name='Freq')
Do_1 = m.FV(value=Do_const, lb=100, ub=200,name='Do_1')
Do_2 = m.FV(value=Do_const, lb=100, ub=200,name='Do_2')
p1 = m.Var(value=1, lb=1, ub=20,name='p1')
p2 = m.Var(value=1, lb=1, ub=20,name='p2')
t1 = m.FV(value=copper_thickness, lb=0.5, ub=2,name='t1')
t2 = m.FV(value=copper_thickness, lb=0.5, ub=2,name='t2')
d = m.Var(value=1, lb=2, ub=5,name='d')
N = m.Var(value=4, lb=2, ub=10, integer=True,name='N')
w1 = m.Var(value=13.87, lb=5, ub=50,name='w1')
w2 = m.Var(value=15.06, lb=5, ub=50,name='w2')
Di_1 = m.Var(value=50, lb=5, ub=50,name='Di_1')
Di_2 = m.Var(value=50, lb=5, ub=50,name='Di_2')


# Define constraints
m.Equation(Di_1 == Di_2)
m.Equation(Di_1 <= Do_const)
m.Equation(Di_2 <= Do_const)
m.Equation(t1 + t2 + d <= sys_thickness)
m.Equation(p1 == (Do_1 - Di_1 - (3 * np.pi * N + 4 * np.pi - 1) / (2 * np.pi) * w1) * 2 * np.pi / (3 * np.pi * N - 1))
m.Equation(p2 == (Do_2 - Di_2 - (3 * np.pi * N + 4 * np.pi - 1) / (2 * np.pi) * w2) * 2 * np.pi / (3 * np.pi * N - 1))
m.Equation(p1 <= w1)
m.Equation(p2 <= w2)
# m.Equation(resonant_freq <= (Freq+0.1*Freq)*1e6)
# m.Equation(resonant_freq >= (Freq-0.1*Freq)*1e6)

# Define all variables for prediction
variables = [Freq, Do_1, Do_2, d, p1, p2, N, w1, w2, Di_1, Di_2, t1, t2]

# Define model from Tensorflow/Keras trained model
pred_Q_factor = Gekko_NN_TF(model, mma, m, n_output=1, activationFxn='relu').predict(variables)

# Define the GEKKO objective function
m.Maximize(pred_Q_factor)


# Variable options
Freq.STATUS = 0
Do_1.STATUS = 0
Do_2.STATUS = 0
t1.STATUS = 0
t2.STATUS = 0

m.options.SOLVER = 1  # APOPT in MINLP
m.options.IMODE = 2  # Steady-state parameter estimation mode
m.solve(debug=0, disp=True)
m.open_folder()  # use this to look at infeasibilities.txt
    
# Print the results
params_max_Q = []
print("\nOptimized design:")
for i in range(len(variables)):
    print(f'\t{features[i]}: {variables[i].value[0]: .3f}')
    params_max_Q.append(variables[i].value)
params_max_Q = np.array(params_max_Q).reshape(1, len(variables))
print(f"Maximized Q_factor: {-m.options.OBJFCNVAL}")
print(f"Maximized Q_factor: {model.predict(params_max_Q)[0][0]: .3f}")
print('Gekko Solvetime:', m.options.SOLVETIME, 'sec')
import time
min_time = time.strftime("%H:%M:%S", time.gmtime(m.options.SOLVETIME))
print('Gekko Solvetime:', min_time)


## SOLUTION ##
 ----------------------------------------------------------------
 APMonitor, Version 1.0.3
 APMonitor Optimization Suite
 ----------------------------------------------------------------
 
 
 --------- APM Model Size ------------
 Each time step contains
   Objects      :           97
   Constants    :            0
   Variables    :          309
   Intermediates:           98
   Connections  :          291
   Equations    :          301
   Residuals    :          203
 
 Number of state variables:            498
 Number of total equations: -          396
 Number of slack variables: -            5
 ---------------------------------------
 Degrees of freedom       :             97
 
 ----------------------------------------------
 Model Parameter Estimation with APOPT Solver
 ----------------------------------------------
Iter:     1 I:  0 Tm:      0.43 NLPi:   34 Dpth:    0 Lvs:    0 Obj: -1.72E+05 Gap:  0.00E+00
 Successful solution
 
 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :   0.440399999999499      sec
 Objective      :   -172228.826336312     
 Successful solution
 ---------------------------------------------------
 


Optimized design:
    Freq [MHz]:  13.560
    Do_1 [mm]:  200.000
    Do_2 [mm]:  200.000
    Gap_Distance [mm]:  5.000
    Pitch_1 [mm]:  3.039
    Pitch_2 [mm]:  4.433
    Turns:  10.000
    Width_1 [mm]:  6.229
    Width_2 [mm]:  5.000
    Di_1 [mm]:  50.000
    Di_2 [mm]:  50.000
    Thickness_1 [mm]:  0.500
    Thickness_2 [mm]:  0.500
Maximized Q_factor: 172228.82634
WARNING:tensorflow:6 out of the last 6 calls to <function TensorFlowTrainer.make_predict_function.<locals>.one_step_on_data_distributed at 0x175d93370> triggered tf.function retracing. Tracing is expensive and the excessive number of tracings could be due to (1) creating @tf.function repeatedly in a loop, (2) passing tensors with different shapes, (3) passing Python objects instead of tensors. For (1), please define your @tf.function outside of the loop. For (2), @tf.function has reduce_retracing=True option that can avoid unnecessary retracing. For (3), please refer to https://www.tensorflow.org/guide/function#controlling_retracing and https://www.tensorflow.org/api_docs/python/tf/function for  more details.
[1m1/1[0m [32m━━━━━━━━━━━━━━━━━━━━[0m[37m[0m [1m0s[0m 41ms/step
Maximized Q_factor:  783.136
Gekko Solvetime: 0.4404 sec
Gekko Solvetime: 00:00:00

2. 现在，我想添加一个更严格的约束：

添加约束

但是，gekko无法用这个约束解决这个问题。我想知道如何调试它。如果您知道任何信息，请告诉我。

代码如下所示： 

## 2nd CODE ##
import numpy as np
import pandas as pd
from gekko import GEKKO
from gekko.ML import Gekko_NN_TF, CustomMinMaxGekkoScaler
from tensorflow.keras.models import load_model
import warnings
warnings.filterwarnings("ignore")

# Q-factor prediction from trained TF model
model_path = 'model_path'
model = load_model(model_path)

# CSV file reading
data_path = 'datapath'
data = pd.read_csv(data_path)

# Define output feature of the model
target = ['Q_factor']
# Basically, the input features are [Freq, Do_1, Do_2, d, p1, p2, N, w1, w2, Di_1, Di_2, t1, t2]
features = [col for col in data.columns if col not in target]

# Find max/min of scaler for Gekko_NN_TF parameters
s = CustomMinMaxGekkoScaler(data, features, target)
mma = s.minMaxValues()

# Initialize GEKKO
m = GEKKO(remote=False)

# Design Constants
epsilon = m.Const(value=8.854e-12, name='epsilon')
epsilon_0 = m.Const(value=1.000, name='epsilon_0')
mu = m.Const(value=1.2566e-6, name='mu')

# Design Requirements
freq_const = 13.56
Do_const = 200
sys_thickness = 200
copper_thickness = 0.5

# Define variables with bounds
Freq = m.FV(value=freq_const, lb=1, ub=30,name='Freq')
Do_1 = m.FV(value=Do_const, lb=100, ub=200,name='Do_1')
Do_2 = m.FV(value=Do_const, lb=100, ub=200,name='Do_2')
p1 = m.Var(value=7.06, lb=1, ub=30,name='p1')
p2 = m.Var(value=5.46, lb=1, ub=30,name='p2')
t1 = m.FV(value=copper_thickness, lb=0.5, ub=2,name='t1')
t2 = m.FV(value=copper_thickness, lb=0.5, ub=2,name='t2')
d = m.Var(value=4.38, lb=2, ub=10,name='d')
N = m.Var(value=4, lb=2, ub=10, integer=True,name='N')
w1 = m.Var(value=13.87, lb=5, ub=50,name='w1')
w2 = m.Var(value=15.06, lb=5, ub=50,name='w2')
Di_1 = m.Var(value=50, lb=5, ub=100,name='Di_1')
Di_2 = m.Var(value=50, lb=5, ub=100,name='Di_2')


# Average of some variables for C_coil and L_coil equation
Do = m.Intermediate((Do_1+Do_2)/2, name='Do')
Di = m.Intermediate((Di_1+Di_2)/2, name='Di')
w = m.Intermediate((w1 + w2) / 2, name='avg_w')
t = m.Intermediate((t1 + t2) / 2, name='avg_t')

# Define intermediates variables (reduce model complexity)
C_coil = m.Intermediate(epsilon * epsilon_0 * np.pi * w * 1e-3 * N * (Do + Di) * 1e-3 / (2 * d * 1e-3) * (1 + d / (np.pi * w) * m.log(2 * d / (np.pi * w)) 
                                                                                                          + d / (np.pi * w) * m.log(1 + 2 * d / (np.pi * w) + 2 * m.sqrt((t / d) + (t / d) ** 2))), name='C_coil')
L_coil = m.Intermediate(mu * (N ** 2) * (Do + Di) * 1e-3 / 4 * m.log(2.46 * (Do + Di) / (Do - Di) + 0.2 * ((Do - Di) / (Do + Di)) ** 2), name='L_coil')
resonant_freq = m.Intermediate(1/(2*np.pi*m.sqrt(L_coil*C_coil)), name='resonant_Freq')

# Define constraints
m.Equation(Di_1 == Di_2)
m.Equation(Di_1 <= Do_const)
m.Equation(Di_2 <= Do_const)
m.Equation(t1 + t2 + d <= sys_thickness)
m.Equation(p1 == (Do_1 - Di_1 - (3 * np.pi * N + 4 * np.pi - 1) / (2 * np.pi) * w1) * 2 * np.pi / (3 * np.pi * N - 1))
m.Equation(p2 == (Do_2 - Di_2 - (3 * np.pi * N + 4 * np.pi - 1) / (2 * np.pi) * w2) * 2 * np.pi / (3 * np.pi * N - 1))
m.Equation(p1 <= w1)
m.Equation(p2 <= w2)
## Added constraint
m.Equation(resonant_freq <= (Freq+0.1*Freq)*1e6)
m.Equation(resonant_freq >= (Freq-0.1*Freq)*1e6)

# Define all variables for prediction
variables = [Freq, Do_1, Do_2, d, p1, p2, N, w1, w2, Di_1, Di_2, t1, t2]

# Define model from Tensorflow/Keras trained model
pred_Q_factor = m.Param(Gekko_NN_TF(model, mma, m, n_output=1, activationFxn='relu').predict(variables), lb=0, ub=2000, name='pred_Q_factor')

# Define the GEKKO objective function
m.Maximize(pred_Q_factor)

# Variable options
Freq.STATUS = 0
Do_1.STATUS = 0
Do_2.STATUS = 0
t1.STATUS = 0
t2.STATUS = 0

m.options.SOLVER = 1  # APOPT in MINLP
m.solver_options = [
    'minlp_maximum_iterations 100', 
    'minlp_max_iter_with_int_sol 100', 
    'minlp_as_nlp 0', 
    'nlp_maximum_iterations 50', 
    'minlp_branch_method 1', 
    'minlp_print_level 8',  
    'minlp_integer_tol 1e-4', 
    'minlp_gap_tol 0.001'
]
m.options.IMODE = 2  # Steady-state parameter estimation mode
m.solve(debug=0, disp=True)

m.open_folder()  # use this to look at infeasibilities.txt
    
# Print the results
params_max_Q = []
print("\nOptimized design:")
for i in range(len(variables)):
    print(f'\t{features[i]}: {variables[i].value[0]: .3f}')
    params_max_Q.append(variables[i].value)
params_max_Q = np.array(params_max_Q).reshape(1, len(variables))
print(f"Maximized Q_factor: {-m.options.OBJFCNVAL}")
print(f"Maximized Q_factor: {model.predict(params_max_Q)[0][0]: .3f}")
print('Gekko Solvetime:', m.options.SOLVETIME, 'sec')

**它显示 **Warning: no more possible trial points and no integer solution** 和 @error: Solution Not Found

这个结果是： ** 

 ----------------------------------------------------------------
 APMonitor, Version 1.0.3
 APMonitor Optimization Suite
 ----------------------------------------------------------------
 
 
 --------- APM Model Size ------------
 Each time step contains
   Objects      :           97
   Constants    :            3
   Variables    :          312
   Intermediates:          105
   Connections  :          291
   Equations    :          310
   Residuals    :          205
 
 Number of state variables:            500
 Number of total equations: -          398
 Number of slack variables: -            7
 ---------------------------------------
 Degrees of freedom       :             95
 
 ----------------------------------------------
 Model Parameter Estimation with APOPT Solver
 ----------------------------------------------
Iter:     1 I: -9 Tm:      1.72 NLPi:   51 Dpth:    0 Lvs:    0 Obj:  0.00E+00 Gap:       NaN
 **Warning: no more possible trial points and no integer solution**
 Maximum iterations
 
 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :    1.74070000000211      sec
 Objective      :    97000.0000000000     
 Unsuccessful with error code            0
 ---------------------------------------------------
 
 Creating file: infeasibilities.txt
 @error: Solution Not Found


Optimized design:
    Freq [MHz]:  13.560
    Do_1 [mm]:  200.000
    Do_2 [mm]:  200.000
    Gap_Distance [mm]:  2.000
    Pitch_1 [mm]:  1.000
    Pitch_2 [mm]:  1.000
    Turns:  4.000
    Width_1 [mm]:  13.870
    Width_2 [mm]:  15.060
    Di_1 [mm]:  50.000
    Di_2 [mm]:  50.000
    Thickness_1 [mm]:  0.500
    Thickness_2 [mm]:  0.500
Maximized Q_factor: -97000.0
[1m1/1[0m [32m━━━━━━━━━━━━━━━━━━━━[0m[37m[0m [1m0s[0m 46ms/step
Maximized Q_factor:  902.816
Gekko Solvetime: 1.7407 sec
Gekko Solvetime: 00:00:01

我希望获得最高 pred_Q_factor 和其他参数值。

你能让我知道这里有什么问题吗？以及我该如何调试这个？

问题在于添加的共振频率约束使得问题不可行，再加上 MINLP 优化问题的复杂性。下面是导致警告和错误的原因以及如何调试和修复 Gekko 代码的逐步说明。

1. 约束的不可行性

附加约束： 

m.Equation(resonant_freq <= (Freq+0.1*Freq)*1e6)
m.Equation(resonant_freq >= (Freq-0.1*Freq)*1e6)

对共振频率施加了一个窄带，这可能与其他约束（特别是与 C_coil 、 L_coil 和这些变量的输入参数相关的约束）相冲突。换句话说，可能不存在满足所有约束的解，从而导致优化问题不可行。

2. MINLP 的复杂性

混合整数非线性规划 (MINLP) 问题难以解决，因为解空间是离散的（由于整数变量 N ）和非凸的（由于非线性方程）。这意味着求解器可能难以找到满足所有约束并找到全局最优解的可行点。

调试和修复

A. 分析不可行性

  • 检查 infeasibilities.txt ： 如 Gekko 输出中所述，此文件提供有关不可行约束和可能原因的见解。它将确定哪些方程（或其组合）导致了问题。
  • 放宽约束： 尝试略微放宽共振频率约束，以查看这是否会产生可行解。例如：

python m.Equation(resonant_freq <= (Freq+0.15*Freq)*1e6) m.Equation(resonant_freq >= (Freq-0.15*Freq)*1e6) * 简化模型： 如果可能，请尝试简化计算 C_coil 和 L_coil 的方程。有时，高度非线性表达式会导致数值问题。

B. 提高求解器性能

  • 初始猜测： 为变量提供一个接近预期最优解的良好初始猜测可以帮助求解器更快地找到可行解。尝试根据对系统的了解初始化变量。
  • 缩放： 确保所有变量都具有相似的数量级。缩放可以改善求解器的数值稳定性和性能。
  • 求解器选项： 调整求解器选项（如所做的那样）可以影响其性能。尝试不同的选项，例如：
    • 增加迭代次数： minlp_maximum_iterations 、 nlp_maximum_iterations
    • 调整容差： minlp_integer_tol 、 minlp_gap_tol
    • 尝试不同的分支方法： minlp_branch_method
    • 增加打印级别以获取更多信息： minlp_print_level
  • 不同的求解器： 如果 APOPT 无法找到解，请尝试 Gekko 中可用的其他 MINLP 求解器，例如 IPOPT（对于非线性问题）或尝试全局求解器（如果局部最优解不可接受）。

C. 逐步调试

  • 从更简单的问题开始： 首先尝试在没有整数变量 N 的情况下求解问题（使其成为 NLP）。这可以帮助确定不可行性是来自非线性约束还是来自整数变量。
  • 逐渐添加约束： 逐个添加约束，并在每次添加后求解问题。这将有助于查明哪个特定约束导致了不可行性。

其他提示

  • 检查模型定义： 确保所有方程式、变量和约束都在 Gekko 中正确定义和实现。仔细检查单位和数量级。
  • 使用 m.open_folder() ： 利用此函数检查求解器生成的中间文件，这可以提供有关求解过程的见解。

请记住，解决 MINLP 问题可能具有挑战性，通常需要反复试验和调整求解器设置才能获得良好的解决方案。

标签：python,gekko
From： 78840772

相关文章

  • ArcPro (3.2+) Python 脚本工具中从 .atbx Toolbox 相对导入本地模块
    我设置了一个库和关联的ArcGISToolbox，以便：/root├──Toolbox.atbx├──mylib│└──my_function.py├──my_tools│└──my_gp_script.py我将代码存储库的开发克隆保存在公司共享服务器上的一个位置，并在GitHub上托管一份副本。当我进行更新时，我会......
  • Python vs. R：揭秘机器学习领域的双璧
    一、引言1.1背景介绍随着大数据和人工智能技术的飞速发展，机器学习已经成为了一个热门领域。在机器学习领域，Python和R是两种广泛使用的编程语言。Python因其简洁易读的语法和强大的库支持，成为了最受欢迎的编程语言之一。而R则以其强大的统计分析和数据可视化能力，在统计学......
  • Windows10 安装编译后的 pysqlcipher3-1.2.1 基于 Python 3.8.10
    Windows10安装编译后的pysqlcipher3-1.2.1基于Python3.8.10本文主要是将直接安装编译后的文件，不一定的成功，但是可以尝试使用，若无法直接安装，请参考编译过程，自行编译安装，编译过程见这里安装pysqlcipher3这里用32位举例因为64位安装完全相同，只需要把对应的位数换成64......
  • 【Python】Python基础语法知识点汇集
    Python是一种高级的、解释型的编程语言，以其清晰的语法和代码可读性而闻名。本篇文章将汇集Python编程的基础语法知识点，为初学者提供一个全面的学习指南。......
  • 启动Python 的内置服务器访问本地图片
    要使用Python的内置服务器访问本地图片并正确地显示在浏览器中，你需要将图片文件放在内置服务器的根目录或其子目录中。以下是详细步骤：1.将图片文件复制到服务器根目录：例如，将zheng.jpeg文件复制到一个特定的目录中（例如，你的项目目录）。假设你将图片文件复制到C:\Users\panda......
  • 17：Python数据类型练习题
    #1获取c1，c2相同的元素列表c1=[11,22,33]c2=[22,33,44]foriinc1:ifiinc2:print(i)#2获取c1中有，c2没有的元素列表foriinc1:ifinotinc2:print(i)#3获取c2中有，c1没有的元素列表foriinc2:ifinotinc1:print(i)#4获......
  • (Jmeter新玩法)Python 调 Jmeter执行参数化jmx脚本
    #Python调Jmeter执行参数化jmx脚本importosfromos.pathimportjoinimporttimeimportrefromstringimportTemplatejmeter_Home=r"F:\softtotal\xxx\bin\jmeter.bat"#jmx文件路径currpath=os.path.dirname(os.path.realpath(__file__))#要运行的jmx脚......
  • python爬虫预备知识三-多进程
    python实现多进程的方法：fork、multiprocessing模块创建多进程。os.fork方法os.fork方法只适合于unix/linux系统，不支持windows系统。fork方法调用一次会返回两次，原因在于操作系统将当前进程（父进程）复制出一份进程（子进程），这两个进程几乎完全相同，fork方法分别在父进程和子进程中......
  • PEP 8 – Python 代码风格指南中文版（七）
    编程建议（2） 定义异常时，应该从Exception类继承，而不是从BaseException类继承。直接从BaseException继承的异常通常是那些几乎不应该被捕获的异常。设计异常层次结构时，应该基于捕获异常的代码可能需要进行的区分，而不是基于异常被抛出的位置。目标是通过编程方式回答“出了......
  • Python-记录一次迭代求和
    importitertoolsdefget_result(hope,list_input):""":paramhope:#期望相加所得参数:paramlist_input:#所有数值:return:"""defgenerate_combination(items,length):forcombinationinitertools.co......

赞助商

阅读排行

网站主页关于我们联系我们网站地图

本网站内容转载自其他媒体，侵权联系[admin##ips99.com]。

Copyright © 2020-2023 IPS99 版权所有 IPS99