可以问一下Gekko中问题吗,使用了m.sum对于变量组中的部分变量求和,infeasibilities.txt中求和对应的变量全部为0是什么原因
其中EQ Number 1和2分别对应m.sum(IP_1[IP_deviceset[0,1]-1:IP_deviceset[0,2]])==IP_deviceset[0,3]和m.sum(IP_2[IP_deviceset[1,1]-1:IP_deviceset[1,2]])==IP_deviceset[1,3],对应变量组中存在不为0的值
完整代码为:
from gekko import GEKKO
import numpy as np
import pandas as pd
def model_sol(initial_data,sim_data,SF_deviceset,IP_deviceset,tp,time,deltT,Garget):
#Constants
n=len(time)
#A/C
Ur=0.16*3600
Cr=16500
Rh_max=2.4
Rc_max=4
#price
J_cons=0.6
J_feedin=0.18
J_sf=0.10
J_sd=0.25
J_ip=0.4
m = GEKKO(remote=False)
T_out =sim_data[:,-1]
Q_int=sim_data[:,3]
Pen=sim_data[:,1]
U_nonflex=sim_data[:,2]
#device variable
SF_1,SF_2,SF_3,SF_4,SF_5,SF_6,SF_7=[],[],[],[],[],[],[]
IP_1,IP_2,IP_3,IP_4,IP_5=[],[],[],[],[]
for i in range(n):
SF_1.append(m.Var(value=initial_data[0,i],lb=0, ub=1, integer=True))
SF_2.append(m.Var(value=initial_data[1,i],lb=0, ub=1, integer=True))
SF_3.append(m.Var(value=initial_data[2,i],lb=0, ub=1, integer=True))
SF_4.append(m.Var(value=initial_data[3,i],lb=0, ub=1, integer=True))
SF_5.append(m.Var(value=initial_data[4,i],lb=0, ub=1, integer=True))
SF_6.append(m.Var(value=initial_data[5,i],lb=0, ub=1, integer=True))
SF_7.append(m.Var(value=initial_data[6,i],lb=0, ub=1, integer=True))
IP_1.append(m.Var(value=initial_data[7,i],lb=0, ub=1, integer=True))
IP_2.append(m.Var(value=initial_data[8,i],lb=0, ub=1, integer=True))
IP_3.append(m.Var(value=initial_data[9,i],lb=0, ub=1, integer=True))
IP_4.append(m.Var(value=initial_data[10,i],lb=0, ub=1, integer=True))
IP_5.append(m.Var(value=initial_data[11,i],lb=0, ub=1, integer=True))
# Create auxiliary binary variables to represent the activation of each sum constraint
# 'activation' will be 1 if the sum constraint is activated, 0 otherwise
activation_S1 = [m.Var(lb=0, ub=1, integer=True) for i in range(SF_deviceset[0,2]-SF_deviceset[0,3]-SF_deviceset[0,1]+1)]
activation_S2 = [m.Var(lb=0, ub=1, integer=True) for i in range(SF_deviceset[1,2]-SF_deviceset[1,3]-SF_deviceset[1,1]+1)]
activation_S3 = [m.Var(lb=0, ub=1, integer=True) for i in range(SF_deviceset[2,2]-SF_deviceset[2,3]-SF_deviceset[2,1]+1)]
activation_S4 = [m.Var(lb=0, ub=1, integer=True) for i in range(SF_deviceset[3,2]-SF_deviceset[3,3]-SF_deviceset[3,1]+1)]
activation_S5 = [m.Var(lb=0, ub=1, integer=True) for i in range(SF_deviceset[4,2]-SF_deviceset[4,3]-SF_deviceset[4,1]+1)]
activation_S6 = [m.Var(lb=0, ub=1, integer=True) for i in range(SF_deviceset[5,2]-SF_deviceset[5,3]-SF_deviceset[5,1]+1)]
activation_S7 = [m.Var(lb=0, ub=1, integer=True) for i in range(SF_deviceset[6,2]-SF_deviceset[6,3]-SF_deviceset[6,1]+1)]
Rc = m.Array(m.Var,n,value=0, lb=0, ub=Rc_max,integer=False)
Rh = m.Array(m.Var,n,value=0, lb=0, ub=Rh_max,integer=False)
k_PV=m.Array(m.Var,n,value=1, lb=0, ub=1,integer=False)
#state_variables
P_SF = m.Array(m.Var, n, lb=0,integer=False)
P_IP = m.Array(m.Var, n, lb=0,integer=False)
P_SD = m.Array(m.Var, n, lb=0,integer=False)
S = m.Array(m.Var, n, lb=0,integer=False)
Tr = m.Array(m.Var, n+1,value=20, lb=18,ub=28,integer=False)
G_cons = m.Array(m.Var, n, lb=0,integer=False)
G_feedin = m.Array(m.Var, n, lb=0,integer=False)
#obj variables
s1=m.Array(m.Var, n, lb=0,integer=False)
#Intermediate
#balance
m.Equation(m.sum(IP_1[IP_deviceset[0,1]-1:IP_deviceset[0,2]])==IP_deviceset[0,3])
m.Equation(m.sum(IP_2[IP_deviceset[1,1]-1:IP_deviceset[1,2]])==IP_deviceset[1,3])
m.Equation(m.sum(IP_3[IP_deviceset[2,1]-1:IP_deviceset[2,2]])==IP_deviceset[2,3])
m.Equation(m.sum(IP_4[IP_deviceset[3,1]-1:IP_deviceset[3,2]])==IP_deviceset[3,3])
m.Equation(m.sum(IP_5[IP_deviceset[4,1]-1:IP_deviceset[4,2]])==IP_deviceset[4,3])
for i in range(SF_deviceset[0,1]-1,SF_deviceset[0,2]-SF_deviceset[0,3]):
m.Equation(m.sum(SF_1[i:i+SF_deviceset[0,3]])== SF_deviceset[0,3]*activation_S1[i-SF_deviceset[0,1]+1])
for i in range(SF_deviceset[1,1]-1,SF_deviceset[1,2]-SF_deviceset[1,3]):
m.Equation(m.sum(SF_2[i:i+SF_deviceset[1,3]])== SF_deviceset[1,3]*activation_S2[i-SF_deviceset[1,1]+1])
for i in range(SF_deviceset[2,1]-1,SF_deviceset[2,2]-SF_deviceset[2,3]):
m.Equation(m.sum(SF_3[i:i+SF_deviceset[2,3]])== SF_deviceset[2,3]*activation_S3[i-SF_deviceset[2,1]+1])
for i in range(SF_deviceset[3,1]-1,SF_deviceset[3,2]-SF_deviceset[3,3]):
m.Equation(m.sum(SF_4[i:i+SF_deviceset[3,3]])== SF_deviceset[3,3]*activation_S4[i-SF_deviceset[3,1]+1])
for i in range(SF_deviceset[4,1]-1,SF_deviceset[4,2]-SF_deviceset[4,3]):
m.Equation(m.sum(SF_5[i:i+SF_deviceset[4,3]])== SF_deviceset[4,3]*activation_S5[i-SF_deviceset[4,1]+1])
for i in range(SF_deviceset[5,1]-1,SF_deviceset[5,2]-SF_deviceset[5,3]):
m.Equation(m.sum(SF_6[i:i+SF_deviceset[5,3]])== SF_deviceset[5,3]*activation_S6[i-SF_deviceset[5,1]+1])
for i in range(SF_deviceset[6,1]-1,SF_deviceset[6,2]-SF_deviceset[6,3]):
m.Equation(m.sum(SF_7[i:i+SF_deviceset[6,3]])== SF_deviceset[6,3]*activation_S7[i-SF_deviceset[6,1]+1])
for i in range(n):
m.Equation(Tr[i+1]==(Ur/Cr)*(T_out[i]-Tr[i])*deltT+(Rh[i]-Rc[i]+Q_int[i])*3600*deltT/Cr)
m.Equation(Rh[i]*Rc[i]<=0)
m.Equation(S[i]==-1*k_PV[i]*Pen[i])
m.Equation(P_SF[i]==SF_deviceset[0,4]*SF_1[i]+SF_deviceset[1,4]*SF_2[i]+SF_deviceset[2,4]*SF_3[i]+SF_deviceset[3,4]*SF_4[i]+SF_deviceset[4,4]*SF_5[i]+SF_deviceset[5,4]*SF_6[i]+SF_deviceset[6,4]*SF_7[i])
m.Equation(P_IP[i]==IP_deviceset[0,4]*IP_1[i]+IP_deviceset[1,4]*IP_2[i]+IP_deviceset[2,4]*IP_3[i]+IP_deviceset[3,4]*IP_4[i]+IP_deviceset[4,4]*IP_5[i])
m.Equation(P_SD[i]==Rh[i]+Rc[i])
m.Equation(G_cons[i]-G_feedin[i]==P_SF[i]+P_IP[i]+P_SD[i]+U_nonflex[i]-S[i])
m.Equation(G_cons[i]*G_feedin[i]<=0)
m.Equation(s1[i]==(J_cons*G_cons[i]-J_feedin*G_feedin[i]+J_sf*P_SF[i]+J_ip*P_IP[i]+J_sd*P_SD[i])**2*deltT**2+(tp[i]*(G_cons[i]-G_feedin[i]-Garget[i]))**2)
m.Obj(m.sum(s1))
#options
#m.options.MAX_ITER=100000
#m.options.IMODE=3#
# m.options.NODES=3
#m.solver_options = ['minlp_gap_tol 0.1',\
# 'minlp_maximum_iterations 1000',\
# 'minlp_max_iter_with_int_sol 500',\
# 'minlp_branch_method 1',\
# 'nlp_maximum_iterations 100']
#m.options.COLDSTART = 1
m.options.SOLVER=1#APORT
#m.options.RTOL=1e-6
#m.open_folder() # open folder if remote=False to see infeasibilities.txt
m.solve(disp=True) # solve
Gcons_val=np.zeros(n)
Gfeedin_val=np.zeros(n)
Psf_val=np.zeros(n)
Pip_val=np.zeros(n)
Psd_val=np.zeros(n)
Pucl_val=np.zeros(n)
S_val=np.zeros(n)
for i in range(n):
Gcons_val[i]=G_cons[i].value[0]
Gfeedin_val[i]=G_feedin[i].value[0]
Psf_val[i]=P_SF[i].value[0]
Pip_val[i]=P_IP[i].value[0]
Psd_val[i]=P_SD[i].value[0]
Pucl_val[i]=U_nonflex[i]
S_val[i]=S[i].value[0]
cost_val=float(np.sum(((J_cons*Gcons_val-J_feedin*Gfeedin_val+J_sf*Psf_val+J_ip*Pip_val+J_sd*Psd_val)*deltT)**2))
print(cost_val)
return Psf_val,Pip_val,Psd_val,Pucl_val,S_val,Gcons_val,Gfeedin_val,cost_val