将两个功能与Gekko中的条件if开关结合使用

编程问答 2022-06-12

问题描述

对于低温（300-1000K）和高温（1000-3000K），我有两个热力学关系。如果要在Gekko中同时使用这两种方法，如何将它们组合成一个可以在优化问题中使用的相关性？

这是一段Python代码，可计算300K至3000K的低温或高温关系。

Heat capacity for Methane

import numpy as np
import matplotlib.pyplot as plt

T = np.linspace(300.0,3000.0,50)

a_lo = np.array([ 5.15,-1.37E-02,4.92E-05,-4.85E-08,1.67E-11])
a_hi = np.array([7.49E-02,1.34E-02,-5.73E-06,1.22E-09,-1.02E-13])
i_lo = np.where(np.logical_and(T>=300.0,T<1000.0))
i_hi = np.where(np.logical_and(T>=1000.0,T<=3000.0))
cp = np.zeros(50)
Rg = 8.314 # J/mol-K
cp[i_lo] = a_lo[0] + a_lo[1]*T[i_lo] + a_lo[2]*T[i_lo]**2.0 + \
             a_lo[3]*T[i_lo]**3.0 + a_lo[4]*T[i_lo]**4.0
cp[i_hi] = a_hi[0] + a_hi[1]*T[i_hi] + a_hi[2]*T[i_hi]**2.0 + \
             a_hi[3]*T[i_hi]**3.0 + a_hi[4]*T[i_hi]**4.0
cp *= Rg

plt.plot(T,cp,'k-',lw=5)
plt.plot(T[i_lo],cp[i_lo],'.',color='orange')
plt.plot(T[i_hi],cp[i_hi],color='red')
plt.xlabel('Temperature (K)'); plt.grid()
plt.ylabel(r'$CH_4$ Heat Capacity $\left(\frac{J}{mol-K}\right)$')
plt.show()

我尝试在构建模型时使用条件（if）语句，但它仅使用从初始化值中选择的相关性。如果温度T在我的模型中是一个变量，我希望它根据温度变量切换到另一个。

解决方法

在优化或仿真问题中有几种使用条件函数的方法。第一种方法并不精确，但可以通过使用三次样条在采样点之间创建插值的方法来作为合适的近似值（请参见方法＃1 ）。第二种方法是精确的，但是需要带有if2()的具有互补约束（MPCC）的数学程序或带有if3()的Integer Switch变量（请参见方法＃2 ）。在Logical Conditions in Optimization的“设计优化课程”页面上讨论了这两种方法。

import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO

# CH4 Heat capacity parameters (LO: 300-1000K,HI: 1000K-3000K)
a_lo = np.array([ 5.15,-1.37E-02,4.92E-05,-4.85E-08,1.67E-11])
a_hi = np.array([7.49E-02,1.34E-02,-5.73E-06,1.22E-09,-1.02E-13])
Rg = 8.314 # J/mol-K
m  = GEKKO()

# Approach #1: Cubic Spline
def cp1(T):
    if T>=300 and T<=1000:
        a = a_lo
    elif T>1000 and T<=3000:
        a = a_hi
    else:
        raise Exception('Temperature ' + str(T) + ' out of range')
    cp = (a[0]+a[1]*T+a[2]*T**2.0+a[3]*T**3.0+a[4]*T**4.0)*Rg
    return cp
# Calculate cp at 50 pts
T = np.linspace(300.0,3000.0,50)
cp = [cp1(Ti) for Ti in T]
x1 = m.Var(lb=300,ub=3000); y1 = m.Var()
m.cspline(x1,y1,T,cp)

# Approach #2: Gekko conditional statements
def cp2(a,T):
    return (a[0]+a[1]*T+a[2]*T**2.0+a[3]*T**3.0+a[4]*T**4.0)*Rg
x2  = m.Var(lb=300,ub=3000)
y2a = m.Intermediate(cp2(a_lo,x2)); 
y2b = m.Intermediate(cp2(a_hi,x2)); 
y2  = m.if3(x2-1000,y2a,y2b)

m.Equation(y1==80)
m.Equation(y2==80)

m.solve()

print('Find Temperature where cp=80 J/mol-K')
print(x1.value[0],x2.value[0])
conditional-statementsgekkomathematical-optimizationpython