GEKKO - 如何修复 Python Gekko Max Equation 错误 - 元素数

Posted

技术标签:

【中文标题】GEKKO - 如何修复 Python Gekko Max Equation 错误 - 元素数【英文标题】:GEKKO - How to fix Python Gekko Max Equation error - number of elements 【发布时间】:2021-11-20 21:45:22 【问题描述】:

我使用 Gekko 优化函数开发了一个脚本。下面的脚本针对许多元素运行。我测试了 20 和 47 个细胞(shapefile 数据集)的优化算法,脚本实现了解决方案。但是,当我测试更大的数据集时,例如包含 160 个元素,会显示以下错误消息:

“APM 模型错误:字符串 > 15000 个字符

考虑将线分解成多个方程”

我阅读了一些解决此问题的建议。我尝试使用 m.sum,但问题仍然存在。

拜托,你能帮我解决这个问题吗?

请在我们转移链接下方找到下载包含 47 个单元格和 160 个单元格的数据集。

https://wetransfer.com/downloads/64cc631237adacc926c67f56124b327a20210928212223/d8a2d7

谢谢

亚历山大。

import geopandas as gpd
import time
import csv
from gekko import GEKKO
import numpy as np
import math
import pandas as pd

m = GEKKO()


A = -0.00000536 
B = -0.0000291 
E = 0.4040771 
r = 0.085 


input_path = 'D:/Alexandre/shapes/Threats/Prototype/BHO50k/Velhas_BHO50k1summ4_47cells.shp'


output_folder = 'D:/Alexandre/shapes/Threats/Prototype/Small_area/resultados'


input_layer = gpd.read_file(input_path)

input_layer = input_layer[
    ['cocursodag', 'cobacia', 'nuareacont', 'nudistbact', 'D0c', 'Ki0', 'Kj0', 'nuareamont', 'deltai', 'It',
     'cost_op_BR', 'Ii_ub', 'Itj', 'cj', 'deltaj2']]

input_layer = input_layer.astype('cobacia': 'string', 'cocursodag': 'string')


count_input_feat = input_layer.shape[0]


row=count_input_feat 
col=10 


input_cobacia = 
ubi = 
numareacont = 
Ki0 = 
Kj0 = 
X = 
deltai2 = 
ai = 
aj = 
D0 = 
Itj = 
It = 
deltaj = 


for row1 in input_layer.iterrows():
    i = row1[0]

    input_cobacia[i] = row1[1]['cobacia'] 
    Ki0[i] = row1[1]['Ki0']+0.001 
    Kj0[i] = row1[1]['Kj0'] 
    X[i] = row1[1]['nuareamont'] 
    deltai2[i]  = row1[1]['deltai'] 
    ai[i] = 5423304*(pow(X[i],-0.1406852)) 
    aj[i] = row1[1]['cj']*100 + row1[1]['cost_op_BR']*100  
    ubi[i] = row1[1]['Ii_ub'] 
    numareacont[i] = row1[1]['nuareacont'] 
    D0[i] = row1[1]['D0c'] 
    It[i] = row1[1]['It'] 
    Itj[i] = row1[1]['Itj'] 
    if Itj[i]<1: 
        deltaj[i] = row1[1]['deltaj2'] * 0.0001
    elif Itj[i]<2:
        deltaj[i] = row1[1]['deltaj2'] * 0.0001
    else:
        deltaj[i] = row1[1]['deltaj2'] * 0.0001


Ii = m.Array(m.Var, (row, col))
Ij = m.Array(m.Var, (row, col))


for i in range(row):
    for j in range(col):
        if It[i] == 0:
            Ii[i, j].value = 0
            Ii[i, j].lower = 0
            Ii[i, j].upper = 5
            Ij[i,j].value = 0
            Ij[i,j].lower = 0
            Ij[i,j].upper = numareacont[i]*0.05*Itj[i]/3.704545

        elif It[i] <= 2:
            Ii[i, j].value = 0
            Ii[i, j].lower = 0
            Ii[i, j].upper = 10
            Ij[i, j].value = 0
            Ij[i, j].lower = 0
            Ij[i, j].upper = numareacont[i]*0.05*Itj[i]/3.704545

        elif It[i] <= 2.5:
            Ii[i, j].value = 0
            Ii[i, j].lower = 0
            Ii[i, j].upper = 15
            Ij[i, j].value = 0
            Ij[i, j].lower = 0
            Ij[i, j].upper = numareacont[i]*0.05*Itj[i]/3.704545

        elif It[i] <= 3:
            Ii[i, j].value = 0
            Ii[i, j].lower = 0
            Ii[i, j].upper = 15
            Ij[i, j].value = 0
            Ij[i, j].lower = 0
            Ij[i, j].upper = numareacont[i]*0.05*Itj[i]/3.704545

        else:
            Ii[i,j].value = 0
            Ii[i,j].lower = 0
            Ii[i,j].upper = 20
            Ij[i,j].value = 0
            Ij[i,j].lower = 0
            Ij[i,j].upper = numareacont[i]*0.05*Itj[i]/3.704545


Ki = m.Array(m.Var, (row, col))
Kj = m.Array(m.Var, (row, col))
indicator = m.Array(m.Var, (row, col))
p = 2


numerator = m.Array(m.Var, (row, col))
denominator = m.Array(m.Var, (row, col))
for row2 in input_layer.iterrows():

    input_cobacia2 = row2[1]['cobacia']
    input_cocursodag = row2[1]['cocursodag']
    input_distance = row2[1]['nudistbact']

    numerator = 0
    denominator = 0

 
    exp = f"cobacia > 'input_cobacia2' and cocursodag.str.startswith('input_cocursodag')"

    for j in range(col):
        for target_feat in input_layer.query(exp).iterrows(): 
            i=target_feat[0]
            target_green_area = Ij[i,j]
            target_distance = target_feat[1]['nudistbact']
            distance = float(target_distance) - float(input_distance)

            idw = 1 / (distance + 1) ** p
            numerator += target_green_area * idw
            denominator += idw



        i=row2[0]
        sum = Ij[i,j]

        if denominator:
            indicator[i,j] = numerator / denominator + sum
        else:
            indicator[i,j] = sum


D0F = m.Array(m.Var, (row, col)) 

for i in range(row): 
    def constraintD0(x):
        return x - 0.2
    for j in range(col): 
        if j == 0: 
            m.fix(Ki[i,j],val = Ki0[i])
            Ki[i,j].lower = 0
            Ki[i,j].upper = 500000
            m.fix(Kj[i,j], val = Kj0[i])
            Kj[i,j].lower = 0
            Kj[i,j].upper = 100000
            m.Equation(D0F[i, j] == A * Ki[i, j] + B * Kj[i, j] + E) 
            D0[i] = D0F[i, j]

        else:
            D0F[i,j].lower = 0
            D0F[i, j].upper = 1
            Ki[i,j].lower = 0
            Ki[i,j].upper = 500000
            Kj[i, j].lower = 0
            Kj[i, j].upper = 100000

            m.Equation(Ki[i,j] - Ki[i,j-1] == Ii[i,j] - deltai2[i] * Ki[i,j-1]) 
            m.Equation(Kj[i,j] - Kj[i,j-1] == Ij[i,j] + deltaj[i] * Kj[i,j-1]+indicator[i,j]) 
            m.Equation(D0F[i,j] == A*Ki[i,j] + B*Kj[i,j] + E)
            m.Equation(D0F[i,j]<=D0[i])

dep = 1 / (1+r) 


z1 = m.sum([m.sum([pow(dep, j)*(ai[i]*Ii[i,j]+aj[i]*Ij[i,j]) for i in range(row)]) for j in range(col)])

# Objective
m.Obj(z1)

m.options.IMODE = 3

m.options.SOLVER = 3

m.options.DIAGLEVEL = 1

m.options.REDUCE=3



try:

    m.solve()    # solve

    # Outputs
    output_Ki = pd.DataFrame(columns=['cobacia'] + list(range(col)))
    output_Kj = pd.DataFrame(columns=['cobacia'] + list(range(col)))
    output_Ii = pd.DataFrame(columns=['cobacia'] + list(range(col)))
    output_Ij = pd.DataFrame(columns=['cobacia'] + list(range(col)))
    output_D0 = pd.DataFrame(columns=['cobacia'] + list(range(col)))
    output_ai = pd.DataFrame(columns=['cobacia'] + list(range(col)))
    output_aj = pd.DataFrame(columns=['cobacia'] + list(range(col)))

    for i in range(row):
        for j in range(col):
            print(Ki)


            output_Ii.loc[i, 'cobacia'] = input_cobacia[i]
            output_Ii.loc[i, j] = Ii[i,j].value[0]

            output_Ij.loc[i, 'cobacia'] = input_cobacia[i]
            output_Ij.loc[i, j] = Ij[i,j].value[0]

            output_Ki.loc[i, 'cobacia'] = input_cobacia[i]
            output_Ki.loc[i, j] = Ki[i,j].value[0]

            output_Kj.loc[i, 'cobacia'] = input_cobacia[i]
            output_Kj.loc[i, j] = Kj[i,j].value[0]

            output_D0.loc[i, 'cobacia'] = input_cobacia[i]
            output_D0.loc[i, j] = D0F[i, j].value[0]

            output_ai.loc[i, 'cobacia'] = input_cobacia[i]
            output_ai.loc[i, j] = ai[i]

            output_aj.loc[i, 'cobacia'] = input_cobacia[i]
            output_aj.loc[i, j] = aj[i]

    df_outputIi = pd.DataFrame(output_Ii)
    df_outputIj = pd.DataFrame(output_Ij)
    df_outputKi = pd.DataFrame(output_Ki)
    df_outputKj = pd.DataFrame(output_Kj)
    df_outputD0 = pd.DataFrame(output_D0)
    df_outputai = pd.DataFrame(output_ai)
    df_outputaj = pd.DataFrame(output_aj)


    with pd.ExcelWriter('output.xlsx') as writer:
        df_outputIi.to_excel(writer, sheet_name="resultado Ii")
        df_outputIj.to_excel(writer, sheet_name="resultado Ij")
        df_outputKi.to_excel(writer, sheet_name="resultado Ki")
        df_outputKj.to_excel(writer, sheet_name="resultado Kj")
        df_outputD0.to_excel(writer, sheet_name="resultado D0")
        df_outputai.to_excel(writer, sheet_name="ai")
        df_outputaj.to_excel(writer, sheet_name="aj")

except:
    print('Not successful')
    from gekko.apm import get_file
    print(m._server)
    print(m._model_name)
    f = get_file(m._server,m._model_name,'infeasibilities.txt')
    f = f.decode().replace('\r','')
    with open('infeasibilities.txt', 'w') as fl:
        fl.write(str(f))


for i in range(row):
    for j in range(col):
        print(Ki[i,j].value)
        print(Kj[i,j].value)
        print(D0F[i,j].value)```

【问题讨论】:

【参考方案1】:

通过使用m.open_folder() 打开运行文件夹,查看m.path 中的APMonitor 模型文件gk0_model.apm。从较小的问题开始,随着问题的扩大,观察方程的大小。此限制(每个方程m.Intermediate() 的中间体可用于更有效地表达模型。 Additional information on Intermediates 可在文档中找到。

【讨论】:

以上是关于GEKKO - 如何修复 Python Gekko Max Equation 错误 - 元素数的主要内容,如果未能解决你的问题,请参考以下文章

如何使用 R reticulate 安装 gekko 包?

Gekko(python) 用于单圈时间优化

如何以图形方式验证我的动态优化结果,与 gekko 中的初始条件进行比较

在 Python 中使用 Gekko 求解函数

适用于 Python 的 Gekko 优化套件 - if3 始终 <0

如何将 gekko 对象插入/附加到现有列表/数组?