OMPR - Writing scalable constraints elegently

Question

Just a brief background: There are 10 cement plants producing two kinds of products (OPC, PPC) supplying material to d demand locations by road or by rail (which can be loaded in multiples of x tonnes) whose cost matrix is LC (rows - Source, columns - Demand destinations). Since the source is defined as a combination of plant, product, and transport mode while the constraint is at a plant product level - I came up with this hack to combine constraints accordingly

Below is the working code for the same.

  library(ompr)
  library(ompr.roi)
  library(ROI.plugin.glpk)
  library(listcomp)

  is_connected <- function(i, j) {
    LC[i, j] > 0
  }

  railroutes = grep(":W:",rownames(LC))
  d = demand[demand$DPP %in% colnames(LC),"DEMAND"]
  rakemultiplier = capacity[capacity$SOURCE %in% rownames(LC),"Multiplier"]
  p1OPC = which(grepl(Plants[1],rNames)&grepl("OPC",rNames))
  p2OPC = which(grepl(Plants[2],rNames)&grepl("OPC",rNames))
  p3OPC = which(grepl(Plants[3],rNames)&grepl("OPC",rNames))
  p4OPC = which(grepl(Plants[4],rNames)&grepl("OPC",rNames))
  # p5OPC = which(grepl(Plants[5],rNames)&grepl("OPC",rNames))
  # p6OPC = which(grepl(Plants[6],rNames)&grepl("OPC",rNames))
  p7OPC = which(grepl(Plants[7],rNames)&grepl("OPC",rNames))
  p8OPC = which(grepl(Plants[8],rNames)&grepl("OPC",rNames))
  p9OPC = which(grepl(Plants[9],rNames)&grepl("OPC",rNames))
  p10OPC = which(grepl(Plants[10],rNames)&grepl("OPC",rNames))
  p1PPC = which(grepl(Plants[1],rNames)&grepl("PPC",rNames))
  p2PPC = which(grepl(Plants[2],rNames)&grepl("PPC",rNames))
  p3PPC = which(grepl(Plants[3],rNames)&grepl("PPC",rNames))
  p4PPC = which(grepl(Plants[4],rNames)&grepl("PPC",rNames))
  p5PPC = which(grepl(Plants[5],rNames)&grepl("PPC",rNames))
  p6PPC = which(grepl(Plants[6],rNames)&grepl("PPC",rNames))
  p7PPC = which(grepl(Plants[7],rNames)&grepl("PPC",rNames))
  p8PPC = which(grepl(Plants[8],rNames)&grepl("PPC",rNames))
  p9PPC = which(grepl(Plants[9],rNames)&grepl("PPC",rNames))
  p10PPC = which(grepl(Plants[10],rNames)&grepl("PPC",rNames))
  set.seed(40)
  
  
  a = Sys.time()
  model <- MIPModel() |> 
    add_variable(x[i, j], i = 1:nrow(LC), j = 1:ncol(LC), is_connected(i, j), lb = 0, type = "continuous") |>
    add_variable(y[i], lb = 0, type = "integer", i = railroutes) |>
    add_variable(z1[j],lb = 0,ub = pub,type = "continuous", j = 1:ncol(LC)) |>
    add_variable(z2[j],lb = 0,ub = nub,type = "continuous", j = 1:ncol(LC)) |>
    add_constraint(sum_over(x[i,j],i = 1:nrow(LC),is_connected(i, j))+z1[j]-z2[j] == d[j],j = 1:ncol(LC)) |>
    add_constraint(sum_expr(x[i,j],j = 1:ncol(LC),is_connected(i,j))/rakemultiplier[i] == y[i], i = railroutes) |>
    add_constraint(sum_expr(x[i,j],j = 1:ncol(LC),is_connected(i,j)) == y[i]*rakemultiplier[i], i = railroutes) |>
    add_constraint(sum_over(x[i,j], i = p1OPC, j = 1:ncol(LC),is_connected(i,j)) <= 120000) |>
    add_constraint(sum_over(x[i,j], i = p2OPC, j = 1:ncol(LC),is_connected(i,j)) <= 120000) |>
    add_constraint(sum_over(x[i,j], i = p3OPC, j = 1:ncol(LC),is_connected(i,j)) <= 120000) |>
    add_constraint(sum_over(x[i,j], i = p4OPC, j = 1:ncol(LC),is_connected(i,j)) <= 120000) |>
    # add_constraint(sum_over(x[i,j], i = p5OPC, j = 1:ncol(LC),is_connected(i,j)) <= 0) |>
    # add_constraint(sum_over(x[i,j], i = p6OPC, j = 1:ncol(LC),is_connected(i,j)) <= 0) |>
    add_constraint(sum_over(x[i,j], i = p7OPC, j = 1:ncol(LC),is_connected(i,j)) <= 175000) |>
    add_constraint(sum_over(x[i,j], i = p8OPC, j = 1:ncol(LC),is_connected(i,j)) <= 175000) |>
    add_constraint(sum_over(x[i,j], i = p9OPC, j = 1:ncol(LC),is_connected(i,j)) <= 175000) |>
    add_constraint(sum_over(x[i,j], i = p10OPC, j = 1:ncol(LC),is_connected(i,j)) <= 175000) |>
    add_constraint(sum_over(x[i,j], i = p1PPC, j = 1:ncol(LC),is_connected(i,j)) <= 75000) |>
    add_constraint(sum_over(x[i,j], i = p2PPC, j = 1:ncol(LC),is_connected(i,j)) <= 75000) |>
    add_constraint(sum_over(x[i,j], i = p3PPC, j = 1:ncol(LC),is_connected(i,j)) <= 75000) |>
    add_constraint(sum_over(x[i,j], i = p4PPC, j = 1:ncol(LC),is_connected(i,j)) <= 75000) |>
    add_constraint(sum_over(x[i,j], i = p5PPC, j = 1:ncol(LC),is_connected(i,j)) <= 75000) |>
    add_constraint(sum_over(x[i,j], i = p6PPC, j = 1:ncol(LC),is_connected(i,j)) <= 150000) |>
    add_constraint(sum_over(x[i,j], i = p7PPC, j = 1:ncol(LC),is_connected(i,j)) <= 150000) |>
    add_constraint(sum_over(x[i,j], i = p8PPC, j = 1:ncol(LC),is_connected(i,j)) <= 150000) |>
    add_constraint(sum_over(x[i,j], i = p9PPC, j = 1:ncol(LC),is_connected(i,j)) <= 150000) |>
    add_constraint(sum_over(x[i,j], i = p10PPC, j = 1:ncol(LC),is_connected(i,j)) <= 150000) |>
    set_objective(sum_over(x[i,j]*LC[i,j]+z1[j]+z2[j],i = 1:nrow(LC),j = 1:ncol(LC),is_connected(i, j)) ,"min")
  model
  print(Sys.time() - a)
  a = Sys.time()
  solution1 = model %>% solve_model(with_ROI(solver = "glpk",verbose=T))
  print(Sys.time() - a)
  rm(a)

This works perfectly fine, but one thing I am unhappy about is the fact that this code has to be manually changed when the original data is subsetted. Also when there are some plants that don't produce a certain product, I have to comment on the code - i.e the corresponding rows are numerical(0) (check p5opc,p6opc)

Is there a clever way to combine these expressions into something scalable to x plants taking care of the commenting when the corresponding variable is numerical(0) thus making the code more readable?