Add a new constraint to the model

The complete executable code is as follows:

from gurobipy import Model, GRB, quicksum
import sys
import numpy as np
import time
import os

from DataRead import getdata

from itertools import chain

def MPModel(Data):
    n = Data['n']
    m = Data['m']
    J = Data['J']
    OJ = Data['OJ']
    oj=OJ
    operations_machines = Data['operations_machines']
    operations_times = Data['operations_times']
    largeM = Data['largeM']
    # starttime=[0,1,2,3,4,5,6,7,8]
    D = [[60, 60, 65, 60, 65, 90, 60, 65],
                [65, 65, 85, 75, 60, 75, 85, 80],
                [60, 70, 75, 80, 75, 80, 80, 75],
                [65, 65, 80, 75, 80, 75, 75, 80],
                [60, 65, 65, 60, 75, 80, 80, 75],
                [60, 65, 65, 80, 75, 75, 80, 75],
                [60, 60, 75, 80, 75, 70, 60, 75]]

    # preprocessing
    # obtain the earliest starting time for each operation
    stimeOp = {}
    stimeOpMax = 0
    for i in J:
        for j in OJ[i]:
            if j == 1:
                # stimeOp[i, j] =     starttime[i]
                stimeOp[i, j] = 0
            else:
                opt_time_min = sys.maxsize
                for k in operations_machines[i, j - 1]:
                    if operations_times[i, j - 1, k] < opt_time_min:
                        opt_time_min = operations_times[i, j - 1, k]
                stimeOp[i, j] = stimeOp[i, j - 1] + opt_time_min
            if stimeOpMax < stimeOp[i, j]:
                stimeOpMax = stimeOp[i, j]
    print(stimeOp)
    # the latest completion time for each operation
    ltimeOp = {}
    for i in J:
        for jj in range(len(OJ[i]) - 1, -1, -1):
            j = OJ[i][jj]
            opt_time_min = sys.maxsize
            for k in operations_machines[i, j]:
                if operations_times[i, j, k] < opt_time_min:
                    opt_time_min = operations_times[i, j, k]

            if jj == len(OJ[i]) - 1:
                ltimeOp[i, j] = largeM - opt_time_min
            else:
                ltimeOp[i, j] = ltimeOp[i, j + 1] - opt_time_min
    print(ltimeOp)

    model = Model("model")

    x, y, s = {}, {}, {}
    cmax = model.addVar(lb=stimeOpMax, ub=largeM, vtype="I", name="cmax")

    # define x
    for j in J:  # job
        for i in OJ[j]:  # operation
            # s[j,i]=model.addVar(lb=0,ub=largeM,vtype="I",name="s(%s,%s)"%(j,i))
            s[j, i] = model.addVar(lb=stimeOp[j, i], ub=largeM, vtype="I", name="s(%s,%s)" % (j, i))  # ltimeOp[j,i]
            for k in operations_machines[j, i]:
                x[j, i, k] = model.addVar(lb=0, ub=1, vtype="I", name="x(%s,%s,%s)" % (j, i, k))
    # define y
    for i in J:
        for ip in J:
            for j in OJ[i]:
                for jp in OJ[ip]:
                    y[i, j, ip, jp] = model.addVar(lb=0, ub=1, vtype="I", name="y(%s,%s,%s,%s)" % (i, j, ip, jp))

    # define objective function
    model.setObjective(cmax, GRB.MINIMIZE)
    
    # constraint(3)
    for i in J:
        for j in OJ[i]:
            model.addConstr(sum(x[i, j, k] for k in operations_machines[i, j]) == 1, "assignment(%s,%s)" % (i, j))
    # constraint(4)
    for i in J:
        for j in OJ[i]:
            if j != OJ[i][0]:

                p = None
                current_machine = None
                model.update()
                model.addConstr(s[i, j] >= s[i, j - 1] + quicksum(
                    operations_times[i, j - 1, k] * x[i, j - 1, k] for k in operations_machines[i, j - 1])
                     ,"stime(%s,%s)" % (i, j))
    # constraint(4)  idea1
    # for i in J:
    #     for j in OJ[i]:
    #         if j != OJ[i][0]:

    #             p = None
    #             current_machine = None
    #             model.update()
    #             model.addConstr(s[i, j] >= s[i, j - 1] + quicksum(
    #                 operations_times[i, j - 1, k] * x[i, j - 1, k] for k in operations_machines[i, j - 1])+
    #                 quicksum(D[p*x[i, j - 1, p]][q*x[i, j, q]] for p in operations_machines[i, j - 1] for q in operations_machines[i, j])
    #                 ,"stime(%s,%s)" % (i, j))
    # # constraint(4)  idea2
    # for i in J:
    #     for j in OJ[i]:
    #         if j != OJ[i][0]:

    #             previous_machine = None
    #             current_machine = None
    #             model.update()
    #             for k in operations_machines[i, j-1]:
    #                 if  x[i, j-1, k] > 0.5: # 根据机器分配情况确定前一个操作所在的机器
    #                     previous_machine = k
    #                 break
    #             for k in operations_machines[i, j]:
    #                 if x[i, j, k] > 0.5:  # 根据机器分配情况确定当前操作所在的机器
    #                     current_machine = k
    #                 break
    #             model.addConstr(s[i, j] >= s[i, j - 1] + quicksum(
    #                 operations_times[i, j - 1, k] * x[i, j - 1, k] for k in operations_machines[i, j - 1])+
    #                D[[previous_machine][current_machine]]
    #                ,"stime(%s,%s)" % (i, j))
       # constraint(5) 
    for k in operations_machines[1, 1]    :
        for i in J:
            model.addConstr(sum(x[i, j, k]    for j in OJ[i] ) == 1, "cons_x(%s,%s)" % (i, j))   
    # constraint(6)
    for i in J:
        for ip in J:
            if i < ip:
                for j in OJ[i]:
                    for jp in OJ[ip]:
                        kkp = [k for k in operations_machines[i, j] if k in operations_machines[ip, jp]]
                        if len(kkp):
                            for iipk in kkp:
                                model.addConstr(s[ip, jp] >= s[i, j] + operations_times[i, j, iipk] - largeM * (
                                            3 - y[i, j, ip, jp] - x[i, j, iipk] - x[ip, jp, iipk]),
                                                "cons_6_(%s,%s,%s,%s)" % (i, j, ip, jp))
                                model.addConstr(
                                    s[i, j] >= s[ip, jp] +
                                    operations_times[ip, jp, iipk] - largeM *
                                    (2 + y[i, j, ip, jp] - x[i, j, iipk] - x[ip, jp, iipk]),
                                    "cons_5_(%s,%s,%s,%s)" % (i, j, ip, jp))
    # cmax constraint
    for i in J:
        model.addConstr(cmax >= s[i, OJ[i][len(OJ[i]) - 1]] + quicksum(
            operations_times[i, OJ[i][len(OJ[i]) - 1], k] * x[i, OJ[i][len(OJ[i]) - 1], k] for k in
            operations_machines[i, OJ[i][len(OJ[i]) - 1]]), "cmax_cons_(%s)" % (i))

    model.params.TimeLimit=90
    model.update()
    model.optimize()
    # print(y)

    for j in J:  # job
        for i in OJ[j]:  # operation
            # s[j,i]=model.addVar(lb=0,ub=largeM,vtype="I",name="s(%s,%s)"%(j,i))
           
            for k in operations_machines[j, i]:
                string= ''
                if abs(  x[j, i, k].x-1) <= 0.01:
                    string =j+i+k+s[j,i].x
                    print(j,i,k,s[j,i].x)
                    # print(s[j,i].x)
    # print(s)
    return model
filename='datasss.txt'
Data=getdata(filename)

n=Data['n']
m=Data['m']
J=Data['J']
OJ=Data['OJ']
print(OJ)
operations_machines=Data['operations_machines']
operations_times=Data['operations_times']
largeM=Data['largeM']

mipmodel=MPModel(Data) 

txt document content is:
4 4 1
4 4 3 54 1 34 4 61 2 2 4 3 54 1 34 4 61 2 2 4 3 54 1 34 4 61 2 2 4 3 54 1 34 4 61 2 2
4 4 4 9 1 15 2 89 3 70 4 4 9 1 15 2 89 3 70 4 4 9 1 15 2 89 3 70 4 4 9 1 15 2 89 3 70
4 4 1 38 2 19 3 28 4 87 4 1 38 2 19 3 28 4 87 4 1 38 2 19 3 28 4 87 4 1 38 2 19 3 28 4 87
4 4 1 95 3 34 2 7 4 29 4 1 95 3 34 2 7 4 29 4 1 95 3 34 2 7 4 29 4 1 95 3 34 2 7 4 29

txt explanation：

-in the first line there are (at least) 2 numbers: the first is the number of jobs and the second the number of machines (the 3rd is not necessary, it is the average number of machines per operation)

-Every row represents one job: the first number is the number of operations of that job, the second number (let's say k>=1) is the number of machines that can process the first operation; then according to k, there are k pairs of numbers (machine,processing time) that specify which are the machines and the processing times; then the data for the second operation and so on...

Example:  6x6 instance, txt：

6   6   1   
6   1   3   1   1   1   3   1   2   6   1   4   7   1   6   3   1   5   6   
6   1   2   8   1   3   5   1   5   10  1   6   10  1   1   10  1   4   4   
6   1   3   5   1   4   4   1   6   8   1   1   9   1   2   1   1   5   7   
6   1   2   5   1   1   5   1   3   5   1   4   3   1   5   8   1   6   9   
6   1   3   9   1   2   3   1   5   5   1   6   4   1   1   3   1   4   1   
6   1   2   3   1   4   3   1   6   9   1   1   10  1   5   4   1   3   1   

first row = 6 jobs and 6 machines 1 machine per operation
second row: job 1 has 6 operations, the first operation can be processed by 1 machine that is machine 3 with processing time 1.

************question：
My code is a model of shop scheduling that takes into account many constraints, and now I want to add a new constraint: the cost of transportation between machines (I've already written the distance cost matrix D).
The existing constraint is that the start time of the current operation s[i, j] must be later than the end time of the previous operation s[i, j-1]+ the time taken by the previous operation
There is a distance transport matrix D, representing the time cost between all machines, and the start time of the current operation s[i, j] must be later than the end time of the previous operation s[i, j-1]+ the time taken by the previous operation + the transport time of the machine on which the two operations are located
In the code, x[j, i, k] indicates that the JTH step of job i operates (machining) on machine k, so i must also consider the order of the two machines K-1 and k (k >1), and k and k+1(k<OJ[i]) when adding the distance cost.
I tried to add two types of code in the comment in constraint 4, and I got the following two types of errors:

TypeError: list indices must be integers or slices, not gurobipy.   LinExpr

result: TypeError: '>' not supported between instances of 'Var' and 'float'

I don't know how to reasonably add this constraint