How to use the pyscipopt.Model function in PySCIPOpt

def flp(I, J, M, d,f, c=None, monolithic=False):
    """flp -- model for the capacitated facility location problem
    Parameters:
        - I: set of customers
        - J: set of facilities
        - d[i]: demand for customer i
        - M[j]: capacity of facility j
        - f[j]: fixed cost for using a facility in point j
        - c[i,j]: unit cost of servicing demand point i from facility j
    Returns a model, ready to be solved.
    """

    master = Model("flp-master")
    # creating the problem
    y = {}
    for j in J:
        y["y(%d)"%j] = master.addVar(vtype="B", name="y(%s)"%j)

    if monolithic:
        x = {}
        demand = {}
        capacity = {}
        for j in J:
            for i in I:
                x[i, j] = master.addVar(vtype="C", name="x(%s,%s)" % (i, j))

        for i in I:
            demand[i] = master.addCons(quicksum(x[i, j] for j in J) >= d[i], "Demand(%s)" % i)

def test_reopt(self):

        m = Model()
        m.enableReoptimization()

        x = m.addVar(name="x", ub=5)
        y = m.addVar(name="y", lb=-2, ub=10)


        m.addCons(2 * x + y >= 8)
        m.setObjective(x + y)
        m.optimize()
        print("x", m.getVal(x))
        print("y", m.getVal(y))
        self.assertEqual(m.getVal(x), 5.0)
        self.assertEqual(m.getVal(y), -2.0)

        m.freeReoptSolve()
        m.addCons(y <= 3)

def flp(I, J, M, d,f, c=None, monolithic=False):
    """flp -- model for the capacitated facility location problem
    Parameters:
        - I: set of customers
        - J: set of facilities
        - d[i]: demand for customer i
        - M[j]: capacity of facility j
        - f[j]: fixed cost for using a facility in point j
        - c[i,j]: unit cost of servicing demand point i from facility j
    Returns a model, ready to be solved.
    """

    master = Model("flp-master")
    # creating the problem
    y = {}
    for j in J:
        y["y(%d)"%j] = master.addVar(vtype="B", name="y(%s)"%j)

    if monolithic:
        x = {}
        demand = {}
        capacity = {}
        for j in J:
            for i in I:
                x[i, j] = master.addVar(vtype="C", name="x(%s,%s)" % (i, j))

        for i in I:
            demand[i] = master.addCons(quicksum(x[i, j] for j in J) >= d[i], "Demand(%s)" % i)

- a0: lower bound on the cycle time (x axis)
        - aK: upper bound on the cycle time (x axis)
        - K: number of linear pieces to use in the approximation
    Returns a model, ready to be solved.
    """

    # construct points for piecewise-linear relation, store in a,b
    a,b = {},{}
    delta = float(aK-a0)/K
    for i in I:
        for k in range(K):
            T = a0 + delta*k
            a[i,k] = T                          # abscissa: cycle time
            b[i,k] = F[i]/T + h[i]*d[i]*T/2.    # ordinate: (convex) cost for this cycle time

    model = Model("multi-item, capacitated EOQ")

    x,c,w_ = {},{},{}
    for i in I:
        x[i] = model.addVar(vtype="C", name="x(%s)"%i)  # cycle time for item i
        c[i] = model.addVar(vtype="C", name="c(%s)"%i)  # total cost for item i
        for k in range(K):
            w_[i,k] = model.addVar(ub=1, vtype="C", name="w(%s,%s)"%(i,k)) #todo ??

    for i in I:
        model.addCons(quicksum(w_[i,k] for k in range(K)) == 1)
        model.addCons(quicksum(a[i,k]*w_[i,k] for k in range(K)) == x[i])
        model.addCons(quicksum(b[i,k]*w_[i,k] for k in range(K)) == c[i])

    model.addCons(quicksum(w[i]*d[i]*x[i] for i in I) <= W)

    model.setObjective(quicksum(c[i] for i in I), "minimize")

def scheduling_linear_ordering(J,p,d,w):
    """
    scheduling_linear_ordering: model for the one machine total weighted tardiness problem

    Model for the one machine total weighted tardiness problem
    using the linear ordering formulation

    Parameters:
        - J: set of jobs
        - p[j]: processing time of job j
        - d[j]: latest non-tardy time for job j
        - w[j]: weighted of job j; the objective is the sum of the weighted completion time

    Returns a model, ready to be solved.
    """
    model = Model("scheduling: linear ordering")

    T,x = {},{} # tardiness variable; x[j,k] =1 if job j precedes job k, =0 otherwise
    for j in J:
        T[j] = model.addVar(vtype="C", name="T(%s)"%(j))
        for k in J:
            if j != k:
                x[j,k] = model.addVar(vtype="B", name="x(%s,%s)"%(j,k))

    for j in J:
        model.addCons(quicksum(p[k]*x[k,j] for k in J if k != j) - T[j] <= d[j]-p[j], "Tardiness(%r)"%(j))

        for k in J:
            if k <= j:
                continue
            model.addCons(x[j,k] + x[k,j] == 1, "Disjunctive(%s,%s)"%(j,k))

def prodmix(I,K,a,p,epsilon,LB):
    """prodmix:  robust production planning using soco
    Parameters:
        I - set of materials
        K - set of components
        a[i][k] -  coef. matrix
        p[i] - price of material i
        LB[k] - amount needed for k
    Returns a model, ready to be solved.
    """

    model = Model("robust product mix")

    x,rhs = {},{}
    for i in I:
        x[i] = model.addVar(vtype="C", name="x(%s)"%i)
    for k in K:
        rhs[k] = model.addVar(vtype="C", name="rhs(%s)"%k)

    model.addCons(quicksum(x[i] for i in I) == 1)
    for k in K:
        model.addCons(rhs[k] == -LB[k]+ quicksum(a[i,k]*x[i] for i in I) )
        model.addCons(quicksum(epsilon*epsilon*x[i]*x[i] for i in I) <= rhs[k]*rhs[k])

    model.setObjective(quicksum(p[i]*x[i] for i in I), "minimize")

    model.data = x,rhs
    return model

def mctransp(I,J,K,c,d,M):
    """mctransp -- model for solving the Multi-commodity Transportation Problem
    Parameters:
        - I: set of customers
        - J: set of facilities
        - K: set of commodities
        - c[i,j,k]: unit transportation cost on arc (i,j) for commodity k
        - d[i][k]: demand for commodity k at node i
        - M[j]: capacity
    Returns a model, ready to be solved.
    """

    model = Model("multi-commodity transportation")

    # Create variables
    x = {}

    for (i,j,k) in c:
        x[i,j,k] = model.addVar(vtype="C", name="x(%s,%s,%s)" % (i,j,k))

    # Demand constraints
    for i in I:
        for k in K:
            model.addCons(sum(x[i,j,k] for j in J if (i,j,k) in x) == d[i,k], "Demand(%s,%s)" % (i,k))

    # Capacity constraints
    for j in J:
        model.addCons(sum(x[i,j,k] for (i,j2,k) in x if j2 == j) <= M[j], "Capacity(%s)" % j)

def eld_complete(U,p_min,p_max,d,brk):
    """eld -- economic load dispatching in electricity generation
    Parameters:
        - U: set of generators (units)
        - p_min[u]: minimum operating power for unit u
        - p_max[u]: maximum operating power for unit u
        - d: demand
        - brk[k]: (x,y) coordinates of breakpoint k, k=0,...,K
    Returns a model, ready to be solved.
    """

    model = Model("Economic load dispatching")

    p,F = {},{}
    for u in U:
        p[u] = model.addVar(lb=p_min[u], ub=p_max[u], name="p(%s)"%u)    # capacity
        F[u] = model.addVar(lb=0,name="fuel(%s)"%u)

    # set fuel costs based on piecewise linear approximation
    for u in U:
        abrk = [X for (X,Y) in brk[u]]
        bbrk = [Y for (X,Y) in brk[u]]

        # convex combination part:
        K = len(brk[u])-1
        z = {}
        for k in range(K+1):
            z[k] = model.addVar(ub=1) # do not name variables for avoiding clash

def eld_another(U,p_min,p_max,d,brk):
    """eld -- economic load dispatching in electricity generation
    Parameters:
        - U: set of generators (units)
        - p_min[u]: minimum operating power for unit u
        - p_max[u]: maximum operating power for unit u
        - d: demand
        - brk[u][k]: (x,y) coordinates of breakpoint k, k=0,...,K for unit u
    Returns a model, ready to be solved.
    """
    model = Model("Economic load dispatching")

    # set objective based on piecewise linear approximation
    p,F,z = {},{},{}
    for u in U:
        abrk = [X for (X,Y) in brk[u]]
        bbrk = [Y for (X,Y) in brk[u]]
        p[u],F[u],z[u] = convex_comb_sos(model,abrk,bbrk)
        p[u].lb = p_min[u]
        p[u].ub = p_max[u]

    # demand satisfaction
    model.addCons(quicksum(p[u] for u in U) == d, "demand")

    # objective
    model.setObjective(quicksum(F[u] for u in U), "minimize")

def scf(n,c):
    """scf: single-commodity flow formulation for the (asymmetric) traveling salesman problem
    Parameters:
        - n: number of nodes
        - c[i,j]: cost for traversing arc (i,j)
    Returns a model, ready to be solved.
    """
    model = Model("atsp - scf")

    x,f = {},{}
    for i in range(1,n+1):
        for j in range(1,n+1):
            if i != j:
                x[i,j] = model.addVar(vtype="B", name="x(%s,%s)"%(i,j))
                if i==1:
                    f[i,j] = model.addVar(lb=0, ub=n-1, vtype="C", name="f(%s,%s)"%(i,j))
                else:
                    f[i,j] = model.addVar(lb=0, ub=n-2, vtype="C", name="f(%s,%s)"%(i,j))

    for i in range(1,n+1):
        model.addCons(quicksum(x[i,j] for j in range(1,n+1) if j != i) == 1, "Out(%s)"%i)
        model.addCons(quicksum(x[j,i] for j in range(1,n+1) if j != i) == 1, "In(%s)"%i)

    model.addCons(quicksum(f[1,j] for j in range(2,n+1)) == n-1, "FlowOut")