英文:
Gekko : How to retrieve the result from an array variable and why it does not solve in local
问题
I have two issues : One is not related to the solver but to how to retrieve my results, the other is related to the solver when I work in local
我有两个问题：一个与求解器无关，而与如何检索我的结果有关，另一个与本地工作时与求解器有关
Also I could not find enough intel in the documentation but i was wondering what are the algorithms behind IMODE = 4 and 7 ?
另外，我在文档中找不到足够的信息，但我想知道IMODE = 4和7背后的算法是什么？
The first issue is that I tried everything to retrieve the 100 values of dist for each time to put on an animation. I did not succeed and I do not understand why...
第一个问题是，我尝试了一切来检索每个时间点的100个dist值以用于动画。但我没有成功，我不明白为什么...
Gives me the following error (I spare you all the iteration) :
给我以下错误（我省略了所有迭代）：
ValueError: could not broadcast input array from shape (2,) into shape (100,)
错误：无法从形状（2，）广播输入数组到形状（100，）
The second issue I have is when I work with : m = GEKKO(remote=False)
我遇到的第二个问题是当我使用：m = GEKKO(remote=False)时
Here the result :
这里是结果：
Fortran runtime error: Out of memory
Fortran运行时错误：内存不足
Error: 'results.json' not found.
错误：找不到'results.json'。
I truly find Gekko amazing and I'm working on mastering it for my work !
我真的觉得Gekko很棒，我正在努力掌握它，以便在工作中使用！感谢你们的辛勤工作和奉献！
英文:
I have two issues : One is not related to the solver but to how to retrieve my results, the other is related to the solver when I work in local
Also I could not find enough intel in the documentation but i was wondering what are the algorithms behind IMODE = 4 and 7 ?
The first issue is that I tried everything to retrieve the 100 values of dist for each time to put on an animation. I did not succeed and I do not understand why...
from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation, rc
# Temps de simulation
tf = 100
# Nombre de tranches d'âge et initialisation nb personne dans une tranche
N = 100
x = np.linspace(0, N, N)
y = 10 + 40 * np.exp(-30 * ((x/100) - 0.5)**2) + np.random.normal(0, 1, size=100)
y = (y - np.min(y)) / (np.max(y) - np.min(y)) * 20 + 10
# Nombre de pas de temps
Ndt = tf*1000
# Création de l'instance du modèle GEKKO
m = GEKKO(remote=True)
# m.open_folder()
m.time = np.linspace(0, tf, Ndt)
# Données d'entrée
birth_rate = 0.001 # taux de natalité
death_rate = 0.0005 # taux de mortalité
growth_rate = 0.01 # taux de croissance
# Variables du modèle
dist = m.Array(m.Var, N, lb=0)
birth = m.Array(m.Var, N, lb=0)
death = m.Array(m.Var, N, lb=0)
growth = m.Array(m.Var, N, lb=0)
# Contraintes initiales
for i in range(N):
dist[i].VALUE = y[i]
# Contraintes sur la variable dist
m.Equation([dist[i].dt() == - death[i] + growth[i-1] for i in range(1,N)])
m.Equation([birth[i] == 0 for i in range(1,N)])
m.Equation([death[i] == death_rate * dist[i]*i for i in range(1,N)])
m.Equation([growth[i] == growth_rate * dist[i] for i in range(1,N)])
m.Equation(dist[0].dt() == birth[0] - death[0])
m.Equation(birth[0] == birth_rate * dist[0])
m.Equation(death[0] == death_rate * dist[0])
m.Equation(growth[0] == 0)
# Résolution du modèle
m.options.IMODE = 4
m.solve(disp=True)
# Créer la figure et les axes
fig, ax = plt.subplots()
ydata = np.zeros((N))
# Fonction qui met à jour l'animation
def update(frame):
# Effacer le dessin précédent
ax.clear()
ydata[:] = np.array(dist[frame].VALUE)
# Tracer la ligne
ax.plot(np.arange(len(ydata)), ydata)
# Définir les limites des axes
ax.set_xlim(0, len(dist))
ax.set_ylim(0, np.amax(ydata) )
# Ajouter les étiquettes des axes
ax.set_xlabel('X Label')
ax.set_ylabel('Y Label')
# Ajouter un titre à la figure
ax.set_title(f'Animation Example Temps: {frame*5}')
# Créer l'animation
ani = animation.FuncAnimation(fig, update, frames=Ndt, interval=10)
# Afficher l'animation
ani.save('distribution.gif')
plt.show()
Gives me the following error (I spare you all the iteration) :
Number of Iterations....: 72
(scaled) (unscaled)
Objective...............: 0.0000000000000000e+00 0.0000000000000000e+00
Dual infeasibility......: 0.0000000000000000e+00 0.0000000000000000e+00
Constraint violation....: 5.4521725359095381e-07 5.4521725359095381e-07
Complementarity.........: 0.0000000000000000e+00 0.0000000000000000e+00
Overall NLP error.......: 5.4521725359095381e-07 5.4521725359095381e-07
Number of objective function evaluations = 107
Number of objective gradient evaluations = 39
Number of equality constraint evaluations = 563
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 73
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 72
Total CPU secs in IPOPT (w/o function evaluations) = 0.146
Total CPU secs in NLP function evaluations = 0.212
EXIT: Optimal Solution Found.
The solution was found.
The final value of the objective function is 0.000000000000000E+000
---------------------------------------------------
Solver : IPOPT (v3.12)
Solution time : 0.397200000006706 sec
Objective : 0.000000000000000E+000
Successful solution
---------------------------------------------------
Traceback (most recent call last):
File "c:\Users\yaj\Downloads\Gekko_test.py", line 81, in <module>
ani.save('distribution.gif')
File "C:\Users\yaj\Anaconda3\lib\site-packages\matplotlib\animation.py", line 1068, in save
anim._init_draw() # Clear the initial frame
File "C:\Users\yaj\Anaconda3\lib\site-packages\matplotlib\animation.py", line 1721, in _init_draw
self._draw_frame(frame_data)
File "C:\Users\yaj\Anaconda3\lib\site-packages\matplotlib\animation.py", line 1743, in _draw_frame
self._drawn_artists = self._func(framedata, *self._args)
File "c:\Users\yaj\Downloads\Gekko_test.py", line 62, in update
ydata[:] = np.array(dist[frame].VALUE)
ValueError: could not broadcast input array from shape (2,) into shape (100,)
The second issue I have is when I work with : m = GEKKO(remote=False)
Here the result :
PS C:\Users\yaj\Downloads\H0504-Kinetic Digital Modelling> conda activate base
PS C:\Users\yaj\Downloads\H0504-Kinetic Digital Modelling> & C:/Users/yaj/Anaconda3/python.exe c:/Users/yaj/Downloads/Gekko_test.py
----------------------------------------------------------------
APMonitor, Version 1.0.0
APMonitor Optimization Suite
----------------------------------------------------------------
--------- APM Model Size ------------
Each time step contains
Objects : 0
Constants : 0
Variables : 400
Intermediates: 0
Connections : 0
Equations : 400
Residuals : 400
Error: At line 1545 of file apm.f90
Traceback: not available, compile with -ftrace=frame or -ftrace=full
Fortran runtime error: Out of memory
Error: 'results.json' not found. Check above for additional error details
Traceback (most recent call last):
File "c:\Users\yaj\Downloads\Gekko_test.py", line 50, in <module>
m.solve(disp=True)
File "C:\Users\yaj\Anaconda3\lib\site-packages\gekko\gekko.py", line 2227, in solve
self.load_JSON()
File "C:\Users\yaj\Anaconda3\lib\site-packages\gekko\gk_post_solve.py", line 13, in load_JSON
f = open(os.path.join(self._path,'options.json'))
FileNotFoundError: [Errno 2] No such file or directory: 'C:\\Users\\yaj\\AppData\\Local\\Temp\\tmpnfynca23gk_model0\\options.json'
I truly find Gekko amazing and I'm working on mastering it for my work !
Thank you for all you work and dedication !
答案1
得分: 1
IMODE=4 是一个同时模拟，在这个模式下，所有 100000 个时间点和 4x100=400 个状态 + 1x100 个微分方程一起求解。问题规模是 100000 x 5 x 100 = 50,000,000 个变量和方程。似乎存在将一个长达 100000 行的时间向量传输到远程服务器的问题，因此它只计算一个时间步。然而，将时间向量传输到本地解决方案没有相同的问题，因此尝试创建大问题并耗尽内存。如果是 Windows 计算机，本地求解器是 32 位的，因此进程可以分配的最大内存为 4 GB。在 Linux 和 MacOS 上，是 64 位的本地求解器。然而，使用 5000 万个变量进行大内存分配并不是一个好主意，因为使用非线性规划求解器（如 IPOPT）解决问题需要一段时间。问题似乎是线性的，因此应该在几次迭代内解决。切换到 IMODE=7 可以逐个时间步骤解决微分方程。然而，计算 100000 个时间步需要一些时间。
保持 IMODE=4，减少模拟的时间步数。尝试从每年计算 1000 个时间步减少到每年只计算一个时间步作为起点。APOPT 求解器在 0.27 秒内找到了 50000 个变量和 50000 个方程的解决方案。
以下是源代码：
from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation, rc
# Temps de simulation
tf = 100
# Nombre de tranches d'âge et initialisation nb personne dans une tranche
N = 100
x = np.linspace(0, N, N)
y = 10 + 40 * np.exp(-30 * ((x/100) - 0.5)**2) + np.random.normal(0, 1, size=100)
y = (y - np.min(y)) / (np.max(y) - np.min(y)) * 20 + 10
# Nombre de pas de temps
Ndt = tf+1
# Création de l'instance du modèle GEKKO
m = GEKKO(remote=False)
# m.open_folder()
m.time = np.linspace(0, tf, Ndt)
# Données d'entrée
birth_rate = 0.001 # taux de natalité
death_rate = 0.0005 # taux de mortalité
growth_rate = 0.01 # taux de croissance
# Variables du modèle
dist = m.Array(m.Var, N, lb=0)
birth = m.Array(m.Var, N, lb=0)
death = m.Array(m.Var, N, lb=0)
growth = m.Array(m.Var, N, lb=0)
# Contraintes initiales
for i in range(N):
dist[i].VALUE = y[i]
# Contraintes sur la variable dist
m.Equations([dist[i].dt() == - death[i] + growth[i-1] for i in range(1,N)])
m.Equations([birth[i] == 0 for i in range(1,N)])
m.Equations([death[i] == death_rate * dist[i]*i for i in range(1,N)])
m.Equations([growth[i] == growth_rate * dist[i] for i in range(1,N)])
m.Equation(dist[0].dt() == birth[0] - death[0])
m.Equation(birth[0] == birth_rate * dist[0])
m.Equation(death[0] == death_rate * dist[0])
m.Equation(growth[0] == 0)
# Résolution du modèle
m.options.IMODE = 4
m.options.SOLVER = 1
m.solve(disp=True)
# Créer la figure et les axes
fig, ax = plt.subplots(figsize=(8,4))
# Fonction qui met à jour l'animation
def update(frame):
# Effacer le dessin précédent
ax.clear()
ydata = dist[frame].VALUE
# Tracer la ligne
ax.plot(np.arange(len(ydata)), ydata)
# Définir les limites des axes
ax.set_xlim(0, len(dist))
ax.set_ylim(0, np.amax(ydata) )
# Ajouter les étiquettes des axes
ax.set_xlabel('X Label')
ax.set_ylabel('Y Label')
# Ajouter un titre à la figure
ax.set_title(f'Animation Example Temps: {frame*5}')
# Créer l'animation
ani = animation.FuncAnimation(fig, update, frames=N, interval=10)
# Afficher l'animation
ani.save('distribution.gif')
plt.show()
这是动画的源代码。
英文:
IMODE=4 is a simultaneous simulation where all 100000 time points and 4x100=400 states + 1x100 differentials are solved together. The problem size is 100000 x 5 x 100 = 50,000,000 variables and equations. It appears that there is a problem with transferring a time vector to the remote server that is 100000 lines long so it only calculates one time step. However, it doesn't have the same problem with transferring the time vector for a local solution so it attempts to create the large problem and runs out of memory. If it is a Windows computer, the local solver is 32-bit so the maximum that a process can allocate is 4 GB. On Linux and MacOS, it is a 64-bit local solver. However, a large memory allocation with 50M variables is not a good idea because it would take a while to solve with a nonlinear programming solver such as IPOPT. The problem appears to be linear so it should solve in a couple iterations. Switching to IMODE=7 solves the differential equations one time step at a time. However, 100000 time steps will take a while to compute.
Stay with IMODE=4 and reduce the number of time steps for simulation. Instead of computing 1000 time steps every year, try just one time step for each year as a starting point. APOPT solver finds the solution to 50000 variables and 50000 equations in 0.27 sec.
Number of state variables: 50000
Number of total equations: - 50000
Number of slack variables: - 0
---------------------------------------
Degrees of freedom : 0
----------------------------------------------
Dynamic Simulation with APOPT Solver
----------------------------------------------
Iter Objective Convergence
0 1.36282E-31 7.67439E-01
1 1.36282E-31 7.67439E-01
Successful solution
---------------------------------------------------
Solver : APOPT (v1.0)
Solution time : 0.270599999988917 sec
Objective : 0.000000000000000E+000
Successful solution
---------------------------------------------------
Here is the animation:
with the source code:
from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation, rc
# Temps de simulation
tf = 100
# Nombre de tranches d'âge et initialisation nb personne dans une tranche
N = 100
x = np.linspace(0, N, N)
y = 10 + 40 * np.exp(-30 * ((x/100) - 0.5)**2) + np.random.normal(0, 1, size=100)
y = (y - np.min(y)) / (np.max(y) - np.min(y)) * 20 + 10
# Nombre de pas de temps
Ndt = tf+1
# Création de l'instance du modèle GEKKO
m = GEKKO(remote=False)
# m.open_folder()
m.time = np.linspace(0, tf, Ndt)
# Données d'entrée
birth_rate = 0.001 # taux de natalité
death_rate = 0.0005 # taux de mortalité
growth_rate = 0.01 # taux de croissance
# Variables du modèle
dist = m.Array(m.Var, N, lb=0)
birth = m.Array(m.Var, N, lb=0)
death = m.Array(m.Var, N, lb=0)
growth = m.Array(m.Var, N, lb=0)
# Contraintes initiales
for i in range(N):
dist[i].VALUE = y[i]
# Contraintes sur la variable dist
m.Equations([dist[i].dt() == - death[i] + growth[i-1] for i in range(1,N)])
m.Equations([birth[i] == 0 for i in range(1,N)])
m.Equations([death[i] == death_rate * dist[i]*i for i in range(1,N)])
m.Equations([growth[i] == growth_rate * dist[i] for i in range(1,N)])
m.Equation(dist[0].dt() == birth[0] - death[0])
m.Equation(birth[0] == birth_rate * dist[0])
m.Equation(death[0] == death_rate * dist[0])
m.Equation(growth[0] == 0)
# Résolution du modèle
m.options.IMODE = 4
m.options.SOLVER = 1
m.solve(disp=True)
# Créer la figure et les axes
fig, ax = plt.subplots(figsize=(8,4))
# Fonction qui met à jour l'animation
def update(frame):
# Effacer le dessin précédent
ax.clear()
ydata = dist[frame].VALUE
# Tracer la ligne
ax.plot(np.arange(len(ydata)), ydata)
# Définir les limites des axes
ax.set_xlim(0, len(dist))
ax.set_ylim(0, np.amax(ydata) )
# Ajouter les étiquettes des axes
ax.set_xlabel('X Label')
ax.set_ylabel('Y Label')
# Ajouter un titre à la figure
ax.set_title(f'Animation Example Temps: {frame*5}')
# Créer l'animation
ani = animation.FuncAnimation(fig, update, frames=N, interval=10)
# Afficher l'animation
ani.save('distribution.gif')
plt.show()
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库，让每个人都能够通过互相帮助和分享经验来进步。
评论