利用氢能作为长期储能媒介，并通过优化布局提高氢能网络的利用效率

【实例简介】利用氢能作为长期储能媒介，并通过优化布局提高氢能网络的利用效率
【实例截图】
【核心代码】# Decision Variables

    x = m.addVars(I, K, lb = 0, vtype=GRB.BINARY, name="x")
    h_tr = m.addVars(I, I, T, S, lb = 0, vtype=GRB.CONTINUOUS, name="h_tr")
    epsilon_r = m.addVars(I, T, S, lb = 0, vtype=GRB.CONTINUOUS, name="epsilon_r")
    epsilon_h = m.addVars(I, T, S, lb = 0, vtype=GRB.CONTINUOUS, name="epsilon_h")
    epsilon_dis = m.addVars(I, T, S, lb = 0, vtype=GRB.CONTINUOUS, name="epsilon_dis")
    epsilon_n = m.addVars(I, T, S, lb = 0, vtype=GRB.CONTINUOUS, name="epsilon_n")
    epsilon_cha = m.addVars(I, T, S, lb = 0, vtype=GRB.CONTINUOUS, name="epsilon_cha")
    epsilon_e2h = m.addVars(I, T, S, lb = 0, vtype=GRB.CONTINUOUS, name="epsilon_e2h")
    theta_h = m.addVars(I, T, S, lb = 0, vtype=GRB.CONTINUOUS, name="theta_h")
    theta_b = m.addVars(I, T, S, lb = 0, vtype=GRB.CONTINUOUS, name="theta_b")
    h_se = m.addVars(I, T, S, lb = 0, vtype=GRB.CONTINUOUS, name="h_se")
    w = m.addVars(I, T, S, lb = 0, vtype=GRB.CONTINUOUS, name="w")
    o = m.addVars(I, T, S, lb = 0, vtype=GRB.CONTINUOUS, name="o")
    eps_n = m.addVars(I, T, S, lb = 0, vtype=GRB.CONTINUOUS, name="eps_n")
    epsilon = m.addVars(I, T, S, lb = 0, vtype=GRB.CONTINUOUS, name="epsilon")
    theta_h = m.addVars(I, range(8761), S, lb = 0, vtype=GRB.CONTINUOUS, name="theta_h")
    theta_b = m.addVars(I, range(8761), S, lb = 0, vtype=GRB.CONTINUOUS, name="theta_b")
    psi = m.addVars(I, lb = 0, ub = 1, vtype=GRB.CONTINUOUS, name="psi")

    # Objective Function
    Capex = sum(x[i, k] * (c_ins[k] / sum(1 / (1 delta) ** tau for tau in range(1, lifetime 1))) for i in I for k in K)
    
    fixedOM = quicksum(c_om[k] * x[i, k] for i in I for k in K)

    OandM = quicksum(xi[s] * (
    quicksum(c_tr * h_tr[i, j, t, s] * d[i, j] for i in I for j in I if j != i for t in T)
    quicksum(c_w * w[i, t, s] for i in I for t in T)
    quicksum(c_pen * eps_n[i, t, s] for i in I for t in T) -
    quicksum(rho_h * h_se[i, t, s] for i in I for t in T) -
    quicksum(rho_o * o[i, t, s] for i in I for t in T)
    ) for s in S)

    neg_psi = quicksum(-psi[i] for i in I)

    m.setObjective(Capex fixedOM OandM, GRB.MINIMIZE)
    m.setObjectiveN(neg_psi, index=1, priority=1, weight=1)

    # Constraints
    m.addConstrs(
                (epsilon_r[i, t, s] epsilon_h[i, t, s] epsilon_dis[i, t, s]
                 epsilon_n[i, t, s] - epsilon_cha[i, t, s] - epsilon_e2h[i, t, s] >= chi[i, t, s] for i in I for t in T for s in S),
                name="energy_balance_constraint"
                )

    m.addConstrs(
                (theta_h[i, t, s] == 
                theta_h[i, t-1, s]
                quicksum(h_tr[j, i, t, s] - h_tr[i, j, t, s] for j in I if j != i)
                epsilon_e2h[i, t, s] / lambda_e2h -
                h_se[i, t, s] -
                epsilon_h[i, t, s] / lambda_h2e for i in I for t in T for s in S),
                name='hydrogen_balance_constraint'
                )

    m.addConstrs(
                (h_se[i, t, s] >= psi[i]*kappa[i, t, s] for i in I for t in T for s in S),name="hydrogen_sold_constraint"
                )