如何从数组变量中检索结果以及为什么它不能在本地解决。

huangapple go评论107阅读模式
英文:

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...

  1. from gekko import GEKKO
  2. import numpy as np
  3. import matplotlib.pyplot as plt
  4. from matplotlib import animation, rc
  5. # Temps de simulation
  6. tf = 100
  7. # Nombre de tranches d'âge et initialisation nb personne dans une tranche
  8. N = 100
  9. x = np.linspace(0, N, N)
  10. y = 10 + 40 * np.exp(-30 * ((x/100) - 0.5)**2) + np.random.normal(0, 1, size=100)
  11. y = (y - np.min(y)) / (np.max(y) - np.min(y)) * 20 + 10
  12. # Nombre de pas de temps
  13. Ndt = tf*1000
  14. # Création de l'instance du modèle GEKKO
  15. m = GEKKO(remote=True)
  16. # m.open_folder()
  17. m.time = np.linspace(0, tf, Ndt)
  18. # Données d'entrée
  19. birth_rate = 0.001 # taux de natalité
  20. death_rate = 0.0005 # taux de mortalité
  21. growth_rate = 0.01 # taux de croissance
  22. # Variables du modèle
  23. dist = m.Array(m.Var, N, lb=0)
  24. birth = m.Array(m.Var, N, lb=0)
  25. death = m.Array(m.Var, N, lb=0)
  26. growth = m.Array(m.Var, N, lb=0)
  27. # Contraintes initiales
  28. for i in range(N):
  29. dist[i].VALUE = y[i]
  30. # Contraintes sur la variable dist
  31. m.Equation([dist[i].dt() == - death[i] + growth[i-1] for i in range(1,N)])
  32. m.Equation([birth[i] == 0 for i in range(1,N)])
  33. m.Equation([death[i] == death_rate * dist[i]*i for i in range(1,N)])
  34. m.Equation([growth[i] == growth_rate * dist[i] for i in range(1,N)])
  35. m.Equation(dist[0].dt() == birth[0] - death[0])
  36. m.Equation(birth[0] == birth_rate * dist[0])
  37. m.Equation(death[0] == death_rate * dist[0])
  38. m.Equation(growth[0] == 0)
  39. # Résolution du modèle
  40. m.options.IMODE = 4
  41. m.solve(disp=True)
  42. # Créer la figure et les axes
  43. fig, ax = plt.subplots()
  44. ydata = np.zeros((N))
  45. # Fonction qui met à jour l'animation
  46. def update(frame):
  47. # Effacer le dessin précédent
  48. ax.clear()
  49. ydata[:] = np.array(dist[frame].VALUE)
  50. # Tracer la ligne
  51. ax.plot(np.arange(len(ydata)), ydata)
  52. # Définir les limites des axes
  53. ax.set_xlim(0, len(dist))
  54. ax.set_ylim(0, np.amax(ydata) )
  55. # Ajouter les étiquettes des axes
  56. ax.set_xlabel('X Label')
  57. ax.set_ylabel('Y Label')
  58. # Ajouter un titre à la figure
  59. ax.set_title(f'Animation Example Temps: {frame*5}')
  60. # Créer l'animation
  61. ani = animation.FuncAnimation(fig, update, frames=Ndt, interval=10)
  62. # Afficher l'animation
  63. ani.save('distribution.gif')
  64. plt.show()

Gives me the following error (I spare you all the iteration) :

  1. Number of Iterations....: 72
  2. (scaled) (unscaled)
  3. Objective...............: 0.0000000000000000e+00 0.0000000000000000e+00
  4. Dual infeasibility......: 0.0000000000000000e+00 0.0000000000000000e+00
  5. Constraint violation....: 5.4521725359095381e-07 5.4521725359095381e-07
  6. Complementarity.........: 0.0000000000000000e+00 0.0000000000000000e+00
  7. Overall NLP error.......: 5.4521725359095381e-07 5.4521725359095381e-07
  8. Number of objective function evaluations = 107
  9. Number of objective gradient evaluations = 39
  10. Number of equality constraint evaluations = 563
  11. Number of inequality constraint evaluations = 0
  12. Number of equality constraint Jacobian evaluations = 73
  13. Number of inequality constraint Jacobian evaluations = 0
  14. Number of Lagrangian Hessian evaluations = 72
  15. Total CPU secs in IPOPT (w/o function evaluations) = 0.146
  16. Total CPU secs in NLP function evaluations = 0.212
  17. EXIT: Optimal Solution Found.
  18. The solution was found.
  19. The final value of the objective function is 0.000000000000000E+000
  20. ---------------------------------------------------
  21. Solver : IPOPT (v3.12)
  22. Solution time : 0.397200000006706 sec
  23. Objective : 0.000000000000000E+000
  24. Successful solution
  25. ---------------------------------------------------
  26. Traceback (most recent call last):
  27. File "c:\Users\yaj\Downloads\Gekko_test.py", line 81, in <module>
  28. ani.save('distribution.gif')
  29. File "C:\Users\yaj\Anaconda3\lib\site-packages\matplotlib\animation.py", line 1068, in save
  30. anim._init_draw() # Clear the initial frame
  31. File "C:\Users\yaj\Anaconda3\lib\site-packages\matplotlib\animation.py", line 1721, in _init_draw
  32. self._draw_frame(frame_data)
  33. File "C:\Users\yaj\Anaconda3\lib\site-packages\matplotlib\animation.py", line 1743, in _draw_frame
  34. self._drawn_artists = self._func(framedata, *self._args)
  35. File "c:\Users\yaj\Downloads\Gekko_test.py", line 62, in update
  36. ydata[:] = np.array(dist[frame].VALUE)
  37. 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 :

  1. PS C:\Users\yaj\Downloads\H0504-Kinetic Digital Modelling> conda activate base
  2. PS C:\Users\yaj\Downloads\H0504-Kinetic Digital Modelling> & C:/Users/yaj/Anaconda3/python.exe c:/Users/yaj/Downloads/Gekko_test.py
  3. ----------------------------------------------------------------
  4. APMonitor, Version 1.0.0
  5. APMonitor Optimization Suite
  6. ----------------------------------------------------------------
  7. --------- APM Model Size ------------
  8. Each time step contains
  9. Objects : 0
  10. Constants : 0
  11. Variables : 400
  12. Intermediates: 0
  13. Connections : 0
  14. Equations : 400
  15. Residuals : 400
  16. Error: At line 1545 of file apm.f90
  17. Traceback: not available, compile with -ftrace=frame or -ftrace=full
  18. Fortran runtime error: Out of memory
  19. Error: 'results.json' not found. Check above for additional error details
  20. Traceback (most recent call last):
  21. File "c:\Users\yaj\Downloads\Gekko_test.py", line 50, in <module>
  22. m.solve(disp=True)
  23. File "C:\Users\yaj\Anaconda3\lib\site-packages\gekko\gekko.py", line 2227, in solve
  24. self.load_JSON()
  25. File "C:\Users\yaj\Anaconda3\lib\site-packages\gekko\gk_post_solve.py", line 13, in load_JSON
  26. f = open(os.path.join(self._path,'options.json'))
  27. 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 个方程的解决方案。

以下是源代码:

  1. from gekko import GEKKO
  2. import numpy as np
  3. import matplotlib.pyplot as plt
  4. from matplotlib import animation, rc
  5. # Temps de simulation
  6. tf = 100
  7. # Nombre de tranches d'âge et initialisation nb personne dans une tranche
  8. N = 100
  9. x = np.linspace(0, N, N)
  10. y = 10 + 40 * np.exp(-30 * ((x/100) - 0.5)**2) + np.random.normal(0, 1, size=100)
  11. y = (y - np.min(y)) / (np.max(y) - np.min(y)) * 20 + 10
  12. # Nombre de pas de temps
  13. Ndt = tf+1
  14. # Création de l'instance du modèle GEKKO
  15. m = GEKKO(remote=False)
  16. # m.open_folder()
  17. m.time = np.linspace(0, tf, Ndt)
  18. # Données d'entrée
  19. birth_rate = 0.001 # taux de natalité
  20. death_rate = 0.0005 # taux de mortalité
  21. growth_rate = 0.01 # taux de croissance
  22. # Variables du modèle
  23. dist = m.Array(m.Var, N, lb=0)
  24. birth = m.Array(m.Var, N, lb=0)
  25. death = m.Array(m.Var, N, lb=0)
  26. growth = m.Array(m.Var, N, lb=0)
  27. # Contraintes initiales
  28. for i in range(N):
  29. dist[i].VALUE = y[i]
  30. # Contraintes sur la variable dist
  31. m.Equations([dist[i].dt() == - death[i] + growth[i-1] for i in range(1,N)])
  32. m.Equations([birth[i] == 0 for i in range(1,N)])
  33. m.Equations([death[i] == death_rate * dist[i]*i for i in range(1,N)])
  34. m.Equations([growth[i] == growth_rate * dist[i] for i in range(1,N)])
  35. m.Equation(dist[0].dt() == birth[0] - death[0])
  36. m.Equation(birth[0] == birth_rate * dist[0])
  37. m.Equation(death[0] == death_rate * dist[0])
  38. m.Equation(growth[0] == 0)
  39. # Résolution du modèle
  40. m.options.IMODE = 4
  41. m.options.SOLVER = 1
  42. m.solve(disp=True)
  43. # Créer la figure et les axes
  44. fig, ax = plt.subplots(figsize=(8,4))
  45. # Fonction qui met à jour l'animation
  46. def update(frame):
  47. # Effacer le dessin précédent
  48. ax.clear()
  49. ydata = dist[frame].VALUE
  50. # Tracer la ligne
  51. ax.plot(np.arange(len(ydata)), ydata)
  52. # Définir les limites des axes
  53. ax.set_xlim(0, len(dist))
  54. ax.set_ylim(0, np.amax(ydata) )
  55. # Ajouter les étiquettes des axes
  56. ax.set_xlabel('X Label')
  57. ax.set_ylabel('Y Label')
  58. # Ajouter un titre à la figure
  59. ax.set_title(f'Animation Example Temps: {frame*5}')
  60. # Créer l'animation
  61. ani = animation.FuncAnimation(fig, update, frames=N, interval=10)
  62. # Afficher l'animation
  63. ani.save('distribution.gif')
  64. 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.

  1. Number of state variables: 50000
  2. Number of total equations: - 50000
  3. Number of slack variables: - 0
  4. ---------------------------------------
  5. Degrees of freedom : 0
  6. ----------------------------------------------
  7. Dynamic Simulation with APOPT Solver
  8. ----------------------------------------------
  9. Iter Objective Convergence
  10. 0 1.36282E-31 7.67439E-01
  11. 1 1.36282E-31 7.67439E-01
  12. Successful solution
  13. ---------------------------------------------------
  14. Solver : APOPT (v1.0)
  15. Solution time : 0.270599999988917 sec
  16. Objective : 0.000000000000000E+000
  17. Successful solution
  18. ---------------------------------------------------

Here is the animation:

如何从数组变量中检索结果以及为什么它不能在本地解决。

with the source code:

  1. from gekko import GEKKO
  2. import numpy as np
  3. import matplotlib.pyplot as plt
  4. from matplotlib import animation, rc
  5. # Temps de simulation
  6. tf = 100
  7. # Nombre de tranches d'âge et initialisation nb personne dans une tranche
  8. N = 100
  9. x = np.linspace(0, N, N)
  10. y = 10 + 40 * np.exp(-30 * ((x/100) - 0.5)**2) + np.random.normal(0, 1, size=100)
  11. y = (y - np.min(y)) / (np.max(y) - np.min(y)) * 20 + 10
  12. # Nombre de pas de temps
  13. Ndt = tf+1
  14. # Création de l'instance du modèle GEKKO
  15. m = GEKKO(remote=False)
  16. # m.open_folder()
  17. m.time = np.linspace(0, tf, Ndt)
  18. # Données d'entrée
  19. birth_rate = 0.001 # taux de natalité
  20. death_rate = 0.0005 # taux de mortalité
  21. growth_rate = 0.01 # taux de croissance
  22. # Variables du modèle
  23. dist = m.Array(m.Var, N, lb=0)
  24. birth = m.Array(m.Var, N, lb=0)
  25. death = m.Array(m.Var, N, lb=0)
  26. growth = m.Array(m.Var, N, lb=0)
  27. # Contraintes initiales
  28. for i in range(N):
  29. dist[i].VALUE = y[i]
  30. # Contraintes sur la variable dist
  31. m.Equations([dist[i].dt() == - death[i] + growth[i-1] for i in range(1,N)])
  32. m.Equations([birth[i] == 0 for i in range(1,N)])
  33. m.Equations([death[i] == death_rate * dist[i]*i for i in range(1,N)])
  34. m.Equations([growth[i] == growth_rate * dist[i] for i in range(1,N)])
  35. m.Equation(dist[0].dt() == birth[0] - death[0])
  36. m.Equation(birth[0] == birth_rate * dist[0])
  37. m.Equation(death[0] == death_rate * dist[0])
  38. m.Equation(growth[0] == 0)
  39. # Résolution du modèle
  40. m.options.IMODE = 4
  41. m.options.SOLVER = 1
  42. m.solve(disp=True)
  43. # Créer la figure et les axes
  44. fig, ax = plt.subplots(figsize=(8,4))
  45. # Fonction qui met à jour l'animation
  46. def update(frame):
  47. # Effacer le dessin précédent
  48. ax.clear()
  49. ydata = dist[frame].VALUE
  50. # Tracer la ligne
  51. ax.plot(np.arange(len(ydata)), ydata)
  52. # Définir les limites des axes
  53. ax.set_xlim(0, len(dist))
  54. ax.set_ylim(0, np.amax(ydata) )
  55. # Ajouter les étiquettes des axes
  56. ax.set_xlabel('X Label')
  57. ax.set_ylabel('Y Label')
  58. # Ajouter un titre à la figure
  59. ax.set_title(f'Animation Example Temps: {frame*5}')
  60. # Créer l'animation
  61. ani = animation.FuncAnimation(fig, update, frames=N, interval=10)
  62. # Afficher l'animation
  63. ani.save('distribution.gif')
  64. plt.show()

huangapple
  • 本文由 发表于 2023年2月23日 23:13:29
  • 转载请务必保留本文链接:https://go.coder-hub.com/75546724.html
匿名

发表评论

匿名网友

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen:

确定