我正在尝试用python中的Gekko解决以下问题。

小时
（1） 但是，我遇到了一些错误APM模型错误：字符串>15000个字符考虑将行拆分为多个等式也可能是由于仅使用换行符CR而不是CR LF（对于Windows）或LF（针对MacOS/Linux）。要解决此问题，请使用适当的换行符STOPPING保存APM文件
当我将N增加到更大的数值（如500或600）时，有可能解决这个问题吗？
（2） 我还有一个关于m.min3函数的问题。如果我定义一个变量v，比如v＝m。Var（value=3）并实现m.min3（v，-50），为什么结果是0，而不是-50？

# Import package

from gekko import GEKKO
import numpy as np

# Define parameters

P_CO = 600.0                    # $/tonCO
beta_CO2 = 1.0                  # no unit
P_CO2 = 39.0                    # $/tonCO2eq
E_ref = 3.1022616               # tonCO2eq/tonCO
E_dir = -1.600570692            # tonCO2eq/tonCO
E_indir_others = 0.3339226804   # tonCO2eq/tonCO
E_indir_elec_cons = 18.46607256 # GJ/tonCO
C1_CAPEX = 285695.0             # no unit
C2_CAPEX = 188.42               # no unit
C1_FOX = 82282.0                # no unit
C2_FOX = 24.094                 # no unit
C1_ROX = 4471.5                 # no unit
C2_ROX = 96.034                 # no unit
C1_UOX = 1983.7                 # no unit
C2_UOX = 249.79                 # no unit
r = 0.08                        # discount rate
N = 100                         # number of scenarios
T = 30                          # total time period
GWP_init = 0.338723235          # 2020 Electricity GWP in EU 27 countries
theta_max = 1600000             # Max capacity

# Function to make GWP_EU matrix (TxN matrix)

def Electricity_GWP(GWP_init, n_years, num_episodes):

    GWP_mean = 0.36258224*np.exp(-0.16395611*np.arange(1, n_years+2)) + 0.03091272
    GWP_mean = GWP_mean.reshape(-1,1)
    GWP_Yearly = np.tile(GWP_mean, num_episodes) 

    noise = np.zeros((n_years+1, num_episodes))
    stdev2050 = GWP_mean[-1] * 0.25 
    stdev = np.arange(0, stdev2050 * (1 + 1/n_years), stdev2050/n_years)

    for i in range(n_years+1):
        noise[i,:] = np.random.normal(0, stdev[i], num_episodes) 

    GWP_forecast = GWP_Yearly + noise 

    return GWP_forecast

GWP_EU = Electricity_GWP(GWP_init, T, N) # (T+1)*N matrix
GWP_EU = GWP_EU[1:,:] # T*N matrix

print(np.shape(GWP_EU))

# Build Gekko model

m = GEKKO(remote=False)
theta = m.Array(m.Var, (N,1), lb=0, ub=theta_max)
demand = np.ones((T,1))
demand[0] = 8031887.589
for k in range(1,11):
    demand[k] = demand[k-1] * 1.026 
for k in range(11,21):
    demand[k] = demand[k-1] * 1.016
for k in range(21,T):
    demand[k] = demand[k-1] * 1.011 
demand = 0.12 * demand
demand = np.tile(demand, N) # T*N matrix

print(np.shape(demand))

obj=0
obj = sum(((1/(1+r))**(t+1))*((P_CO*m.min3(demand[t,s], theta[s])) \
        + beta_CO2*P_CO2*m.min3(demand[t,s], theta[s])*(E_ref-E_dir- 
E_indir_others-E_indir_elec_cons*GWP_EU[t,s]) \
        - (C1_CAPEX + C2_CAPEX*theta[s]) - (C1_FOX + C2_FOX*theta[s]) - (C1_ROX + C2_ROX*m.min3(demand[t,s], theta[s])) - (C1_UOX + C2_UOX*m.min3(demand[t,s], theta[s])))
for s in range(N) for t in range(T))/N

m.Maximize(obj)
m.solve()    

print(theta)