How to simplify this model?

Hello,

I'm a Python beginner and I'm trying to program a Time-Space Network in Gurobi, but I don't know if I can simplify it more.

from gurobipy import*
import numpy as np

n = np.array([("000000","004001","005001","050001","      "),
              ("      ","004002","005002","050002","      "),
              ("      ","004003","005003","050003","      "),
              ("      ","004004","005004","050004","      "),
              ("      ","004005","005005","050005","100100")]) 

m = Model("test01")

x_FA1 = m.addVars([(n[a,b]+n[a,c]+"001")
                  for a in range (1,4) 
                  for b in range (1,4) 
                  for c in range (1,4) 
                  if c != b
                  ], vtype=GRB.BINARY, name='x')

x_FA2 = m.addVars([(n[a,b]+n[a,c]+"002")
                  for a in range (2,3)
                  for b in range (1,3) 
                  for c in range (1,3) 
                  if c != b
                  ], vtype=GRB.BINARY, name='x')

x_SA1 = m.addVars([(n[a,a]+n[a,b]+"001")
                  for a in range (0,1) 
                  for b in range (1,3) 
                  ], vtype=GRB.BINARY, name='x')

x_SA2 = m.addVars([(n[a,a]+n[b,c]+"002")
                  for a in range (0,1) 
                  for b in range (1,2)
                  for c in range (1,3)
                  ], vtype=GRB.BINARY, name='x')

x_EA1 = m.addVars([(n[a,b]+n[a,a]+"001")
                  for a in range (4,5) 
                  for b in range (1,3) 
                  ], vtype=GRB.BINARY, name='x')

x_EA2 = m.addVars([(n[a,b]+n[c,c]+"002")
                   for a in range (3,4) 
                  for b in range (1,3)
                  for c in range (4,5)
                  ], vtype=GRB.BINARY, name='x')

x_HA1 = m.addVars([(n[a,c]+n[b,c]+"001")
                   for c in range (1,4) 
                   for a in range (0,4) 
                   for b in range (1,5)
                   if b == a+1
                   ], vtype=GRB.BINARY, name='x')

x_HA2 = m.addVars([(n[a,c]+n[b,c]+"002")
                  for c in range (1,3) 
                  for a in range (1,3) 
                  for b in range (2,4)
                  if b == a+1
                  ], vtype=GRB.BINARY, name='x')
m.update()

m.setObjective(quicksum(x_FA1) + quicksum(x_FA2),GRB.MINIMIZE)

m.addConstr((quicksum(x_SA1)) == 1 ,'c0') 
m.addConstr((quicksum(x_SA2)) == 1 ,'c1') 

m.addConstr(-(quicksum(x_EA1)) == -1 ,'c2') 
m.addConstr(-(quicksum(x_EA2)) == -1 ,'c3') 


#Especially the restricted expressions c4~c28, I can't think of how to write the restricted expressions more efficiently


m.addConstr((x_HA1['004002004003001'] + x_HA2['004002004003002']) <= 1 ,'c4')
m.addConstr((x_HA1['004003004004001'] + x_HA2['004003004004002']) <= 1 ,'c5')
m.addConstr((x_HA1['005002005003001'] + x_HA2['005002005003002']) <= 1 ,'c6') 
m.addConstr((x_HA1['005003005004001'] + x_HA2['005003005004002']) <= 1 ,'c7')


m.addConstr((x_HA1['004001004002001'] - x_SA1['000000004001001']) == 0,"c8")
m.addConstr((x_HA1['005001005002001'] - x_SA1['000000005001001']) == 0,"c9")

m.addConstr((x_HA1['004002004003001'] + x_FA1['004002005002001']+x_FA1['004002050002001'])-(x_HA1['004001004002001']+x_FA1['005002004002001']+x_FA1['050002004002001'])==0,"c10")
m.addConstr((x_HA1['005002005003001'] + x_FA1['005002004002001']+x_FA1['005002050002001'])-(x_HA1['005001005002001']+x_FA1['004002005002001']+x_FA1['050002005002001'])==0,"c11")
m.addConstr((x_HA1['050002050003001'] + x_FA1['050002004002001']+x_FA1['050002005002001'])-(x_FA1['004002050002001']+x_FA1['005002050002001'])==0,"c12")

m.addConstr((x_HA1['004003004004001'] + x_FA1['004003005003001']+x_FA1['004003050003001'])-(x_HA1['004002004003001']+x_FA1['005003004003001']+x_FA1['050003004003001'])==0,"c13")
m.addConstr((x_HA1['005003005004001'] + x_FA1['005003004003001']+x_FA1['005003050003001'])-(x_HA1['005002005003001']+x_FA1['004003005003001']+x_FA1['050003005003001'])==0,"c14")
m.addConstr((x_HA1['050003050004001'] + x_FA1['050003004003001']+x_FA1['050003005003001'])-(x_HA1['050002050003001']+x_FA1['004003050003001']+x_FA1['005003050003001'])==0,"c15")

m.addConstr((x_HA1['004004004005001'] + x_FA1['004004005004001']+x_FA1['004004050004001'])-(x_HA1['004003004004001']+x_FA1['005004004004001']+x_FA1['050004004004001'])==0,"c16")
m.addConstr((x_HA1['005004005005001'] + x_FA1['005004004004001']+x_FA1['005004050004001'])-(x_HA1['005003005004001']+x_FA1['004004005004001']+x_FA1['050004005004001'])==0,"c17")
m.addConstr((x_FA1['050004004004001'] + x_FA1['050004005004001'])-(x_HA1['050003050004001']+x_FA1['004004050004001']+x_FA1['005004050004001'])==0,"c18")

m.addConstr(x_EA1['004005100100001']-x_HA1['004004004005001']==0,"c19")
m.addConstr(x_EA1['005005100100001']-x_HA1['005004005005001']==0,"c20")


m.addConstr(x_HA2['004002004003002'] - x_SA2['000000004002002']==0,"c21")
m.addConstr(x_HA2['005002005003002'] - x_SA2['000000005002002']==0,"c22")
m.addConstr((x_HA2['004003004004002'] + x_FA2['004003005003002'])-(x_HA2['004002004003002'] + x_FA2['005003004003002'])==0,"c23")
m.addConstr((x_HA2['005003005004002'] + x_FA2['005003004003002'])-(x_HA2['005002005003002'] + x_FA2['004003005003002'])==0,"c24")
m.addConstr(x_EA2['004004100100002']-x_HA2['004003004004002']==0,"c25")
m.addConstr(x_EA2['005004100100002']-x_HA2['005003005004002']==0,"c26")


m.addConstr(x_SA2['000000005002002']+x_HA1['004001004002001']+x_HA1['004002004003001']<=2,"c27")


m.addConstr(x_EA2['005004100100002']+x_HA1['004003004004001']+x_FA1['004004050004001']<=2,"c28")

m.optimize()
print("Obj:",m.objVal)
for v in m.getVars():
    if v.x >= 1:
        print("%s:%d"%(v.varName,v.x))

m.write('test01.lp')

Best regards
Sherrody