Decision Optimization

 View Only
Expand all | Collapse all

How to define the constraints 35 and 37

  • 1.  How to define the constraints 35 and 37

    Posted Mon December 26, 2022 05:06 AM
    Edited by System Fri January 20, 2023 04:42 PM
    Dear All,
    I would like to solve a MILP case by using the solver and python from these mathematical constraints and a maximization objective.
    image.png
    However, the result of y ( μ or miu) is 0. (It should be a continuous number more than 0).
    When I checked the constraints one by one in the python code, it seems there are problems with constraint #35 and constraint #38. Perhaps, I defined them incorrectly.
    Technically, since Total_T = 6 then k = [0,1,2,3,4,5] #follow python index
    # constraint 34 is for k = [0,1,2,3,4], namely:
    S_hat[0] - S_hat[1]  ≤ - (v+w)*y
    S_hat[1] - S_hat[2]  ≤ - (v+w)*y
    S_hat[2] - S_hat[3]  ≤ - (v+w)*y
    S_hat[3] - S_hat[4]  ≤ - (v+w)*y
    S_hat[4] - S_hat[5]  ≤ - (v+w)*y
    In addition, # constraint 35 is for k = [5], which is:
    S_hat[5] - (S_hat[1] +1)  ≤ - (v+w)*y
    Meanwhile, # constraint 36 and # constraint 37 are for matching the S_hat and x_in_hat by using the big number (bigM).
    Here is the code.
    import cplex
    from docplex.mp.model import Model
    import numpy as np

    mdl = Model(name='Scheduling')
    inf = cplex.infinity
    bigM= 1000000

    Total_P = 8 # number of places
    Total_T = 6 # number of transitions

    r1 = 10
    r2 = 200
    v = 1
    w = 1 
    h = np.array([v+w,w+r1,v+w,w+r2,v+w,0,0,0]).reshape(Total_P,1)


    AT =np.array([[1,-1,0,0,0,0],
                  [0,1,-1,0,0,0],
                  [0,0,1,-1,0,0],
                  [0,0,0,1,-1,0],
                  [0,0,0,0,1,-1],
                  [0,-1,1,0,0,0],
                  [0,0,0,-1,1,0],
                  [-1,1,-1,1,-1,1]])
    AT_temp = AT.transpose()

    places = np.array(['p1','p2','p3','p4','p5','p6','p7','p8'])
    P_conflict = np.empty((1,Total_P), dtype = object)
    P_zero = np.empty((1,Total_P), dtype = object)

    # define the place without conflict place

    CP = np.count_nonzero(AT, axis=1, keepdims=True) # calculate the nonzero elements for each row

    P_conflict = []
    P_zero = []

    for a in range(0, len(CP)):
        if CP[a].item(0) > 2:
            P_conflict.append(places[a])
        else:
            P_zero.append(places[a])

    y = [ ] # miu

    y = mdl.continuous_var(lb = 0, ub=inf, name='miu') 

    x = np.array(mdl.integer_var_list(Total_T, 0, inf, name='x')).reshape(Total_T,)
    ind_x = np.where(AT[0] == 1)

    def get_index_value(input):
      data = []
      for l in range(len(P_zero)):
        ind_x = np.where(AT[l] == input)
        get_value = x[ind_x]
        data.append(get_value)
      return data


    x_in_hat = get_index_value(1)
    x_out_hat = get_index_value(-1)

    M0 = np.empty((Total_P,1), dtype = object)
    for l in range(Total_P):
        M0[l][0]= mdl.binary_var(name='M0' + str(l+1) + str(',') + str(l+1))

    # Constraint 38 (initial marking constraint

    mdl.add_constraint(M0[1][0] + M0[5][0] == 1)
    mdl.add_constraint(M0[3][0] + M0[6][0] == 1)
    mdl.add_constraint(M0[0][0] + M0[2][0] + M0[4][0] + M0[7][0] == 1)

    M0_temp = M0.reshape(1,Total_P)
    M0_final = M0_temp.tolist()*Total_T

    # constraint 30

    for     l in range(len(P_zero)):
        mdl.add_constraint(x_out_hat[l][0]-x_in_hat[l][0] >= h[l][0]*y - M0[l][0])



    Z = np.empty((Total_T, Total_T), dtype = object)
    for k in range(Total_T):
        for i in range(Total_T):
        #     if k == 0 and i == 0:
        #         Z[0][0]=1
        #     else:
            Z[k][i] = mdl.binary_var(name='Z' + str(k+1) + str(',') + str(i+1))

    storage_ZAT = []
    for k in range(Total_T):
        ZAT = np.matmul(Z[k].reshape(1,Total_T),AT_temp)
        storage_ZAT.append(ZAT)   
    ZAT_final = np.asarray(storage_ZAT).reshape(Total_T,Total_P)

    M = np.empty((Total_T, Total_P), dtype = object)
    for k in range(0,Total_T):
        for l in range (0,Total_P):
            if k == Total_T-1:
                M[Total_T-1][l] = M0_final[0][l]
            else:
                M[k][l] = mdl.integer_var(name='M' + str(k + 1) + str(',') + str(l + 1))

    M_prev = np.empty((Total_T, Total_P), dtype = object)
    if M is not None:
        for k in range(0,Total_T):
            for l in range (0,Total_P):
                if k is not 0:
                    M_prev[k][l] = M[k-1][l]
                else:
                    M_prev[0][l] = M0_final[0][l]

    # constraint 31

    for k in range(Total_T):
      for l in range(Total_P):
            mdl.add_constraint(M[k][l] == M_prev[k][l] + ZAT_final[k][l])

    # constraint 32

    mdl.add_constraints(mdl.sum(Z[k][i] for k in range(Total_T)) == 1 for i in range(Total_T))

    # Constraint 33
    mdl.add_constraints(mdl.sum(Z[k][i] for i in range(Total_T)) == 1 for k in range(Total_T))

    # # # # Parameters
    VW_temp = [[v + w]]
    VW_temp = VW_temp*Total_T
    VW = np.array(VW_temp) 
    S_hat = np.array(mdl.integer_var_list(Total_T, 0, inf, name='S_hat')).reshape(Total_T,1)

    # constraint 34

    for k in range(0,Total_T-1):
         mdl.add_constraint(S_hat[k][0] - S_hat[k+1][0] <= -VW[k][0]*y)

    # Constraint 35

    mdl.add_constraint(S_hat[Total_T-1][0] - (S_hat[0][0] + 1) <= -VW[Total_T-1][0]*y)

    x_temp_hat = x.reshape(Total_T,1)

    # Constraint 36
    for k in range(Total_T):
       for i in range(Total_T):
            mdl.add_constraint(S_hat[k][0] - x_in_hat[i][0] <= (1-Z[k][i])*bigM)

    # Constraint 37
    for k in range(Total_T):
        for i in range(Total_T):
            mdl.add_constraint(S_hat[k][0] - x_in_hat[i][0] >= (Z[k][i]-1)*bigM)



    mdl.maximize(y)
    mdl.print_information()
    solver = mdl.solve()

    if solver is not None:
        mdl.print_solution()
    else:
        print("Solver is error")

    obj_lambda = 1/y.solution_value
    print('obj_lamba=', obj_lambda)

    if mdl.solve():
      print('time', mdl.solve_details.time)

    print('shape VW', np.shape(VW) )
    print('shape S_hat', np.shape(S_hat) )
    Could anyone please help me with what might be incorrect with the python code?
    Thank you.
    Best regards,
    Nicholas
    ​​
    #DecisionOptimization


  • 2.  RE: How to define the constraints 35 and 37

    Posted Wed December 28, 2022 05:37 AM
    Dear Nicholas,

    Given your data, the model is behaving correctly  as if we print values involved in the first constraint #30, we can see:

    x_out_hat[0][0]-x_in_hat[0][0]) >= (h[0][0]*y - M0[0][0]
x_out_hat[0][0] = 64.0
x_in_hat[0][0] = 63.000000000242096
x_out_hat[0][0]-x_in_hat[0][0] = 0.9999999997579039
h[0][0] = 2
M0[0][0] = 0
x_out_hat[1][0]-x_in_hat[1][0]) >= (h[1][0]*y - M0[1][0]
x_out_hat[1][0] = 64.0
x_in_hat[1][0] = 64.0
x_out_hat[1][0]-x_in_hat[1][0] = 0.0
h[1][0] = 11
M0[1][0] = 0
x_out_hat[2][0]-x_in_hat[2][0]) >= (h[2][0]*y - M0[2][0]
x_out_hat[2][0] = 64.00000000000232
x_in_hat[2][0] = 64.0
x_out_hat[2][0]-x_in_hat[2][0] = 2.3163693185779266e-12
h[2][0] = 2
M0[2][0] = 0
x_out_hat[3][0]-x_in_hat[3][0]) >= (h[3][0]*y - M0[3][0]
x_out_hat[3][0] = 64.00000000000232
x_in_hat[3][0] = 64.00000000000232
x_out_hat[3][0]-x_in_hat[3][0] = 0.0
h[3][0] = 201
M0[3][0] = 1.0
x_out_hat[4][0]-x_in_hat[4][0]) >= (h[4][0]*y - M0[4][0]
x_out_hat[4][0] = 65.0
x_in_hat[4][0] = 64.00000000000232
x_out_hat[4][0]-x_in_hat[4][0] = 0.9999999999976836
h[4][0] = 2
M0[4][0] = 0
x_out_hat[5][0]-x_in_hat[5][0]) >= (h[5][0]*y - M0[5][0]
x_out_hat[5][0] = 64.0
x_in_hat[5][0] = 64.0
x_out_hat[5][0]-x_in_hat[5][0] = 0.0
h[5][0] = 0
M0[5][0] = 1.0
x_out_hat[6][0]-x_in_hat[6][0]) >= (h[6][0]*y - M0[6][0]
x_out_hat[6][0] = 64.00000000000232
x_in_hat[6][0] = 64.00000000000232
x_out_hat[6][0]-x_in_hat[6][0] = 0.0
h[6][0] = 0
M0[6][0] = 0

    As x_out_hat is equal to x_in_hat for many indices, it seems the only way for having 0 >= (h * y - M0) for all the values of h and M0  is to have y == 0.
    I hope this helps.
    Cheers




    ------------------------------
    Renaud Dumeur
    ------------------------------

    Original Message