Python:用于求解代数微分方程的GEKKO的快速替代品

Question

为了模拟电网中的潮流，我使用GEKKO来求解我的代数微分方程系统。

对于较小的模拟，它工作得很好，但是连续运行几次，例如在训练强化学习代理时，需要相当长的时间。

谁能建议一个开销更小的求解器，从而加快模拟速度?

我们的GitHub repo链接中有一个系统的小示例:

https://github.com/upb-lea/openmodelica-microgrid-gym/blob/feature_50_SinglePhaseModel/experiments/swing_equation/gekko_freq_volt.py

Answer 1

还有其他的选择，例如assimulo和CasADi。但我还没有找到接近GEKKO的东西。我的建议是，用参数而不是变量来定义模型。这肯定会减少重新初始化模型所需的时间。

Answer 2

在一些基准测试中，Gekko比pyomo快5倍左右。我们还与目前正在使用CasADi的团队合作，通过Gekko加快优化速度。如果它是一种代数建模语言，那么Gekko是最快的选择之一。这里有一些关于Gekko的建议。切换到IMODE=4与m.options.TIME_SHIFT=0，从先前的解决方案作为热启动。这可以显著提高解决方案的速度。另一种选择是通过多线程并行化强化学习代理的函数调用。

import numpy as np
import threading
import time, random
from gekko import GEKKO
class ThreadClass(threading.Thread):
def __init__(self, id, server, ai, bi):
s = self
s.id = id
s.server = server
s.m = GEKKO(remote=False)
s.a = ai
s.b = bi
s.objective = float('NaN')
# initialize variables
s.m.x1 = s.m.Var(1,lb=1,ub=5)
s.m.x2 = s.m.Var(5,lb=1,ub=5)
s.m.x3 = s.m.Var(5,lb=1,ub=5)
s.m.x4 = s.m.Var(1,lb=1,ub=5)
# Equations
s.m.Equation(s.m.x1*s.m.x2*s.m.x3*s.m.x4>=s.a)
s.m.Equation(s.m.x1**2+s.m.x2**2+s.m.x3**2+s.m.x4**2==s.b)
# Objective
s.m.Minimize(s.m.x1*s.m.x4*(s.m.x1+s.m.x2+s.m.x3)+s.m.x3)
# Set global options
s.m.options.IMODE = 3 # steady state optimization
s.m.options.SOLVER = 1 # APOPT solver
threading.Thread.__init__(s)
def run(self):
# Don't overload server by executing all scripts at once
sleep_time = random.random()
time.sleep(sleep_time)
print('Running application ' + str(self.id) + 'n')
# Solve
self.m.solve(disp=False)
# Retrieve objective if successful
if (self.m.options.APPSTATUS==1):
self.objective = self.m.options.objfcnval
else:
self.objective = float('NaN')
self.m.cleanup()
# Select server
server = 'https://byu.apmonitor.com'
# Optimize at mesh points
x = np.arange(20.0, 30.0, 2.0)
y = np.arange(30.0, 50.0, 2.0)
a, b = np.meshgrid(x, y)
# Array of threads
threads = []
# Calculate objective at all meshgrid points
# Load applications
id = 0
for i in range(a.shape[0]):
for j in range(b.shape[1]):
# Create new thread
threads.append(ThreadClass(id, server, a[i,j], b[i,j]))
# Increment ID
id += 1
# Run applications simultaneously as multiple threads
# Max number of threads to run at once
max_threads = 8
for t in threads:
while (threading.activeCount()>max_threads):
# check for additional threads every 0.01 sec
time.sleep(0.01)
# start the thread
t.start()
# Check for completion
mt = 3.0 # max time
it = 0.0 # incrementing time
st = 1.0 # sleep time
while (threading.activeCount()>=1):
time.sleep(st)
it = it + st
print('Active Threads: ' + str(threading.activeCount()))
# Terminate after max time
if (it>=mt):
break
# Wait for all threads to complete
#for t in threads:
#    t.join()
#print('Threads complete')
# Initialize array for objective
obj = np.empty_like(a)
# Retrieve objective results
id = 0
for i in range(a.shape[0]):
for j in range(b.shape[1]):
obj[i,j] = threads[id].objective
id += 1
# plot 3D figure of results
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import cm
import numpy as np
fig = plt.figure()
ax = fig.gca(projection='3d')
surf = ax.plot_surface(a, b, obj, 
rstride=1, cstride=1, cmap=cm.coolwarm, 
vmin = 12, vmax = 22, linewidth=0, antialiased=False)
ax.set_xlabel('a')
ax.set_ylabel('b')
ax.set_zlabel('obj')
ax.set_title('Multi-Threaded GEKKO')
plt.show()

如果有一个DAE积分器，如Sundials中的IDA，那么这也可能是模拟的替代方案。Gekko解决了更高的DAE索引问题，而大多数积分器(如DASSL、DASPK和IDA)仅限于索引1的DAE和可选的索引2的Hessenberg形式。如果可能的话，这通常会给将公式重新排列成这种形式带来很大的负担。