Loop through lpSolve function-r

Question

I have a lpSolve model with 5 periods and I need to repeat the code for each period. Is there a way to automate this process? for instance using a loop (for, while) and conditions? The only thing to consider is that every period output is the next period input for constraints. I currently have 5 periods but it can be 100 periods or more. any help is appreciated.

This is the code (currently with 2 periods for the sake of brevity):

    #Pre
I10 <- as.numeric(0)
I20 <- as.numeric(0)
I30 <- as.numeric(0)
I40 <- as.numeric(0)
I50 <- as.numeric(0)

B10 <- as.numeric(0)
B20 <- as.numeric(0)
B30 <- as.numeric(0)
B40 <- as.numeric(0)
B50 <- as.numeric(0)

#Basic Information (at time 0)
Item <- c(1,2,3,4,5)
Definition <- c("End-item","End-item", "Component","Component","Component")
LeadTime <- c(0,0,0,0,0)
LotSize <- c(replicate(5,"LFL"))
SafetyStock <- c(0,0,0,0,0)
CoefUse_EI1 <- c(0,0,1,1,0)
CoefUse_EI2 <- c(0,0,0,1,1)
BegInventory <- c(I10, I20, I30, I40, I50)
BegBacklog <- c(B10, B20, B30, B40, B50)
ScheduledReceipt <- c(0,0,0,0,0)
BasicInfor <- data.frame(Item,Definition,LeadTime,LotSize,SafetyStock,
                         CoefUse_EI1,CoefUse_EI2,BegInventory, 
                         BegBacklog, ScheduledReceipt)
BasicInfor


#MPS/ DEMNAD FORECAST: (NOTE: component demand = aij *rj)
Period <- c("t0", "t1", "t2", "t3", "t4", "t5")

?rnorm()
D10 <- round(rnorm(1,mean = 10,2));D10
D11 <- round(rnorm(1,mean = 10,2))
D12 <- round(rnorm(1,mean = 10,2))
D13 <- round(rnorm(1,mean = 10,2))
D14 <- round(rnorm(1,mean = 10,2))
D15 <- round(rnorm(1,mean = 10,2))

DemandForecast1 <- c(D10, D11, D12, D13, D14, D15)

D20 <- round(rnorm(1,mean = 10,2))
D21 <- round(rnorm(1,mean = 10,2))
D22 <- round(rnorm(1,mean = 10,2))
D23 <- round(rnorm(1,mean = 10,2))
D24 <- round(rnorm(1,mean = 10,2))
D25 <- round(rnorm(1,mean = 10,2))

DemandForecast2 <- c(D20, D21, D22, D23, D24, D25)

DemandForecast <- rbind(DemandForecast1, DemandForecast2)
colnames(DemandForecast) <- Period
rownames(DemandForecast) <- c("EndItem1", "EndItem2")

DemandForecast

D30 <- CoefUse_EI1[3]* ScheduledReceipt[1] + CoefUse_EI2[3]*ScheduledReceipt[2] ; D30
D40 <- CoefUse_EI1[4]* ScheduledReceipt[1] + CoefUse_EI2[4]*ScheduledReceipt[2] ; D40
D50 <- CoefUse_EI1[5]* ScheduledReceipt[1] + CoefUse_EI2[5]*ScheduledReceipt[2] ; D30

#MODEL

install.packages("lpSolve")
library(lpSolve)

#PERIOD 1
f.obj <- c(3,3,1,1,1,   #I 
           30,30,0,0,0, #B
           0,0,0,0,0
           )

f.con<- matrix(c(
  1,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,
  0,1,0,0,0,0,-1,0,0,0,0,0,0,0,0,
  0,0,1,0,0,0,0,-1,0,0,0,0,0,0,0,
  0,0,0,1,0,0,0,0,-1,0,0,0,0,0,0,
  0,0,0,0,1,0,0,0,0,-1,0,0,0,0,0,
  
  0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,
  0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,
  0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,
  0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
  0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,
  
  0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,
  0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,
  0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,
  0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,
  0,0,0,0,0,0,0,0,0,0,0,0,0,0,1), nrow = 15, byrow = TRUE)

Constraint1.11 <- I10 - B10 - D10 + ScheduledReceipt[1];Constraint1.11
Constraint1.21 <- I20 - B20 - D20 + ScheduledReceipt[2];Constraint1.21
Constraint1.31 <- I30 - B30 - D30 + ScheduledReceipt[3];Constraint1.31
Constraint1.41 <- I40 - B40 - D40 + ScheduledReceipt[4];Constraint1.41
Constraint1.51 <- I50 - B50 - D50 + ScheduledReceipt[5];Constraint1.51

Constraint2.11 <- B10 + D10;Constraint2.11
Constraint2.21 <- B20 + D20;Constraint2.21
Constraint2.31 <- B30 
Constraint2.41 <- B40 
Constraint2.51 <- B50 

Constraint3.11 <- 0
Constraint3.21 <- 0
Constraint3.31 <- 0
Constraint3.41 <- 0
Constraint3.51 <- 0

f.rhs <- c(Constraint1.11, Constraint1.21, Constraint1.31, Constraint1.41,
           Constraint1.51, Constraint2.11, Constraint2.21, Constraint2.31,
           Constraint2.41, Constraint2.51, Constraint3.11, Constraint3.21,
           Constraint3.31, Constraint3.41, Constraint3.51)

f.dir <- c(replicate(5,"="), replicate(5,"<="), replicate(5, ">="))

optimum <- lp(direction = "min", f.obj,f.con,f.dir,f.rhs)

optimum

round(optimum$solution)
optimum$objval

I11 <- optimum$solution[1]; I11
I21 <- optimum$solution[2]
I31 <- optimum$solution[3]
I41 <- optimum$solution[4]
I51 <- optimum$solution[5]
B11 <- optimum$solution[6]
B21 <- optimum$solution[7]
B31 <- optimum$solution[8]
B41 <- optimum$solution[9]
B51 <- optimum$solution[10]
r11 <- optimum$solution[11]
r21 <- optimum$solution[12]
r31 <- optimum$solution[13]
r41 <- optimum$solution[14]
r51 <- optimum$solution[15]

#PERIOD 2: PROBLEM INFEASIBLE SOLUTION
Constraint1.12 <- I11 - B11 - D11 + r11;Constraint1.12
Constraint1.22 <- I21 - B21 - D21 + r21;Constraint1.22
Constraint1.32 <- I31 - B31 - CoefUse_EI1[3]* r11 + r31;Constraint1.32
Constraint1.42 <- I41 - B41 - CoefUse_EI1[4]* r11 - CoefUse_EI2[4]*r21 + r41;Constraint1.42
Constraint1.52 <- I51 - B51 - CoefUse_EI2[5]* r21+ r51;Constraint1.52

Constraint2.12 <- B11 + D11;Constraint2.12
Constraint2.22 <- B21 + D21;Constraint2.22
Constraint2.32 <- B31; Constraint2.32 
Constraint2.42 <- B41; Constraint2.42
Constraint2.52 <- B51; Constraint2.52


f2.rhs <- c(Constraint1.12, Constraint1.22, Constraint1.32, Constraint1.42,
            Constraint1.52, Constraint2.12, Constraint2.22, Constraint2.32,
            Constraint2.42, Constraint2.52, Constraint3.11, Constraint3.21,
            Constraint3.31, Constraint3.41, Constraint3.51)

f2.rhs
length(f2.rhs)


optimum2 <- lp(direction = "min", objective.in =  f.obj, const.mat =  f.con,
              const.dir =  f.dir,const.rhs =  f2.rhs)

optimum2

round(optimum2$solution)
optimum2$objval

I12 <- optimum2$solution[1]; I12
I22 <- optimum2$solution[2]
I32 <- optimum2$solution[3]
I42 <- optimum2$solution[4]
I52 <- optimum2$solution[5]
B12 <- optimum2$solution[6]
B22 <- optimum2$solution[7]
B32 <- optimum2$solution[8]
B42 <- optimum2$solution[9]
B52 <- optimum2$solution[10]
r12 <- optimum2$solution[11]
r22 <- optimum2$solution[12]
r32 <- optimum2$solution[13]
r42 <- optimum2$solution[14]
r52 <- optimum2$solution[15]