? ?????? ?? 2023-07-23 09:49 采纳率: 0%
浏览 38
已结题

Gekko中m.sum求和问题

可以问一下Gekko中问题吗,使用了m.sum对于变量组中的部分变量求和,infeasibilities.txt中求和对应的变量全部为0是什么原因

img

其中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的值

img


完整代码为:

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

  • 写回答

10条回答 默认 最新

  • CSDN专家-sinJack 2023-07-26 15:21
    关注
    获得4.20元问题酬金

    加点日志看看,可能某块代码逻辑有问题

    评论

报告相同问题?

问题事件

  • 系统已结题 7月31日
  • 修改了问题 7月23日
  • 赞助了问题酬金20元 7月23日
  • 创建了问题 7月23日

悬赏问题

  • ¥15 数据量少可以用MK趋势分析吗
  • ¥15 使用VH6501干扰RTR位,CANoe上显示的错误帧不足32个就进入bus off快慢恢复,为什么?
  • ¥15 大智慧怎么编写一个选股程序
  • ¥100 python 调用 cgps 命令获取 实时位置信息
  • ¥15 两台交换机分别是trunk接口和access接口为何无法通信,通信过程是如何?
  • ¥15 C语言使用vscode编码错误
  • ¥15 用KSV5转成本时,如何不生成那笔中间凭证
  • ¥20 ensp怎么配置让PC1和PC2通讯上
  • ¥50 有没有适合匹配类似图中的运动规律的图像处理算法
  • ¥15 dnat基础问题,本机发出,别人返回的包,不能命中