Function constrained optimization in r-Inventory management

Question

I am trying to optimize a function with two constraints (which are also defined as functions) but I keep receiving errors. The main issue is that constraints are long and if I nest them in the objective function, it won't work. Any idea on how to fix it would be appreciated. The current issue is the optim function does not recognize I11 in my constraint1 function

ObjectiveFunction <- function(I){
  
  I11 = I[11]
  I12 = I[12]
  I13 = I[13]
  I14 = I[14]
  I15 = I[15]
  I21 = I[21]
  I22 = I[22]
  I23 = I[23]
  I24 = I[24]
  I25 = I[25]
  I31 = I[31]
  I32 = I[32]
  I33 = I[33]
  I34 = I[34]
  I35 = I[35]
  I41 = I[41]
  I42 = I[42]
  I43 = I[43]
  I44 = I[44]
  I45 = I[45]
  I51 = I[51]
  I52 = I[52]
  I53 = I[53]
  I54 = I[54]
  I55 = I[55]
 
  I61 = I[61]
  I62 = I[62]
  I63 = I[63]
  I64 = I[64]
  I65 = I[65]
  I71 = I[71]
  I72 = I[72]
  I73 = I[73]
  I74 = I[74]
  I75 = I[75]
  I81 = I[81]
  I82 = I[82]
  I83 = I[83]
  I84 = I[84]
  I85 = I[85]
  I91 = I[91]
  I92 = I[92]
  I93 = I[93]
  I94 = I[94]
  I95 = I[95]
  I101 = I[101]
  I102 = I[102]
  I103 = I[103]
  I104 = I[104]
  I105 = I[105] 

  #Objective Function
  Y =
    HoldingCost[1] * (I11 + I12 + I13 + I14 + I15) +
    HoldingCost[2] * (I21 + I22 + I23 + I24 + I25) +
    HoldingCost[3] * (I31 + I32 + I33 + I34 + I35) +
    HoldingCost[4] * (I41 + I42 + I43 + I44 + I45) +
    HoldingCost[5] * (I51 + I52 + I53 + I54 + I55) +
    BacklogCost[1] * (I61 + I62 + I63 + I64 + I65) +
    BacklogCost[2] * (I71 + I72 + I73 + I74 + I75) +
    BacklogCost[3] * (I81 + I82 + I83 + I84 + I85) +
    BacklogCost[4] * (I91 + I92 + I93 + I94 + I95) +
    BacklogCost[5] * (I101 +I102 + I103 + I104 + I105)
  
  if (Constraint1() & Constraint2()) {
    return(Y)
    return(NA)
    
  }
}


Constraint1 <- function(){
  
  I11 - I61 = I10 - B10 - D10 + r1_1  
  I21 - I71 = I20 - B20 - D20 + r2_1 
  I12 - I62 = I11 - I61 - D11 + r1_0 
  I22 - I72 = I21 - I71 - D21 + r2_0 
  
  I13 - I63 = I12 - I62 - D12 + r11 
  I23 - I73 = I22 - I72 - D22 + r21 
      
  I14 - I64 = I13 - I63 - D13 + r12 
  I24 - I74 = I23 - I73 - D23 + r22 
  
  I15 - I65 = I14 - I64 - D14 + r13 
  I25 - I75 = I24 - I74 - D24 + r23 
    
  I31 - I81 = I30 - B30 - D30 + r3_3 
  I41 - I91 = I40 - B40 - D40 + r4_5 
  I51 - I101 = I50 - B50 - D50 + r5_1 
  
  I32 - I82 = I31 - I81 - D31 + r3_2  
  I42 - I92 = I41 - I91 - D41 + r4_4    
  I52 - I102 = I51 - I101 - D51 + r5_0 
        
  I33 - I83 = I32 - I82 - D32 + r3_1 
  I43 - I93 = I42 - I92 - D42 + r4_3 
  I53 - I103 = I52 - I102 - D52 + r51 
    
  I34 - I84 = I33 - I83 - D33 + r3_0 
  I44 - I94 = I43 - I93 - D43 + r4_2 
  I54 - I104 = I53 - I103 - D53 + r52 
      
  I35 - I85 = I34 - I84 - D34 + r31 
  I45 - I95 = I44 - I94 - D44 + r4_1
  I55 - I105 = I54 - I104 - D54 + r53;
}

Constraint2 <- function(){
  
  #Constraint 2
  I61 - B10 <= D10
  I71 - B20 <= D20
  I81 - B30 <= 0
  I91 - B40 <= 0
  I101 - B50 <= 0
  
  I62 - I61 <= D11
  I72 - I71 <= D21
  I82 - I81 <= 0
  I92 - I91 <= 0
  I102 - I101 <= 0
  
  I63 - I62 <= D12
  I73 - I72 <= D22
  I83 - I82 <= 0
  I93 - I92 <= 0
  I103 - I102 <= 0
  
  I64 - I63 <= D13
  I74 - I73 <= D23
  I84 - I83 <= 0
  I94 - I93 <= 0
  I104 - I103 <= 0
  
  I65 - I64 <= D14
  I75 - I74 <= D24
  I85 - I84 <= 0
  I95 - I94 <= 0
  I105 - I104 <= 0
  
}

ParBase <- rep(0,50); ParBase

optim(par = ParBase, ObjectiveFunction)

---

Edit

I am trying to optimize a multi-period inventory control model using a dynamic programming approach and lpSolve. I have split the model into 5 periods. For the first 2 periods, everything works perfectly fine because "r" values are known. However, for period 3 onward, the code cannot work because "I", "B", and "r" are unknown.

I do not know how to solve this issue. I was thinking about the following since there is a relation among "I", "B", and "r" values. If the value of "r" in each period is positive, it equals "I" in that period. else, it is equalized to "B". I was thinking about a loop function but I do not know how to embed such information into the code and if it helps.

This is the code:

    #Basic Information (at time 0)
Item <- c(1,2,3,4,5)
Definition <- c("End-item","End-item", "Component","Component","Component")
LeadTime <- c(1,1,3,5,1)
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(0,0,0,0,0)
BegBacklog <- c(0,0,0,0,0)
ScheduledReceipt <- c(0,0,0,0,0)
ScheduledReceiptT_1 <- c(0,0,0,0,0)
ScheduledReceiptT_2 <- c(0,0,0,0,0)
ScheduledReceiptT_3 <- c(0,0,0,0,0)
ScheduledReceiptT_4 <- c(0,0,0,0,0)
ScheduledReceiptT_5 <- c(0,0,0,0,0)

BasicInfor <- data.frame(Item,Definition,LeadTime,LotSize,SafetyStock,
                         CoefUse_EI1,CoefUse_EI2,BegInventory, 
                         BegBacklog, ScheduledReceiptT_5,
                         ScheduledReceiptT_4,ScheduledReceiptT_3,
                         ScheduledReceiptT_2, ScheduledReceiptT_1,
                         ScheduledReceipt
                         )
BasicInfor

# Predetermined Information
#PARAMETERS : DEMNAD
D10 <- as.numeric(0)
D11 <- as.numeric(20)
D12 <- as.numeric(0)
D13 <- as.numeric(390)
D14 <- as.numeric(0)
D15 <- as.numeric(120)

D20 <- as.numeric(0)
D21 <- as.numeric(20)
D22 <- as.numeric(0)
D23 <- as.numeric(390)
D24 <- as.numeric(0)
D25 <- as.numeric(120)

D30 <- CoefUse_EI1[3]* D10 + CoefUse_EI2[3]*D20 ; D30
D31 <- CoefUse_EI1[3]* D11 + CoefUse_EI2[3]*D21 ; D31
D32 <- CoefUse_EI1[3]* D12 + CoefUse_EI2[3]*D22 ; D32
D33 <- CoefUse_EI1[3]* D13 + CoefUse_EI2[3]*D23 ; D33
D34 <- CoefUse_EI1[3]* D14 + CoefUse_EI2[3]*D24 ; D34
D35 <- CoefUse_EI1[3]* D15 + CoefUse_EI2[3]*D25 ; D35

D40 <- CoefUse_EI1[4]* D10 + CoefUse_EI2[4]*D20 ; D40
D41 <- CoefUse_EI1[4]* D11 + CoefUse_EI2[4]*D21 ; D41
D42 <- CoefUse_EI1[4]* D12 + CoefUse_EI2[4]*D22 ; D42
D43 <- CoefUse_EI1[4]* D13 + CoefUse_EI2[4]*D23 ; D43
D44 <- CoefUse_EI1[4]* D14 + CoefUse_EI2[4]*D24 ; D44
D45 <- CoefUse_EI1[4]* D15 + CoefUse_EI2[4]*D25 ; D45

D50 <- CoefUse_EI1[5]* D10 + CoefUse_EI2[5]*D20 ; D50
D51 <- CoefUse_EI1[5]* D11 + CoefUse_EI2[5]*D21 ; D51
D52 <- CoefUse_EI1[5]* D12 + CoefUse_EI2[5]*D22 ; D52
D53 <- CoefUse_EI1[5]* D13 + CoefUse_EI2[5]*D23 ; D53
D54 <- CoefUse_EI1[5]* D14 + CoefUse_EI2[5]*D24 ; D54
D55 <- CoefUse_EI1[5]* D15 + CoefUse_EI2[5]*D25 ; D55

#SCHEDULED RECEIPT

r1_5 <- ScheduledReceiptT_5[1]
r1_4 <- ScheduledReceiptT_4[1]
r1_3 <- ScheduledReceiptT_3[1]
r1_2 <- ScheduledReceiptT_2[1]
r1_1 <- ScheduledReceiptT_1[1]
r1_0 <- ScheduledReceipt[1]

r2_5 <- ScheduledReceiptT_5[2]
r2_4 <- ScheduledReceiptT_4[2]
r2_3 <- ScheduledReceiptT_3[2]
r2_2 <- ScheduledReceiptT_2[2]
r2_1 <- ScheduledReceiptT_1[2]
r2_0 <- ScheduledReceipt[2]

r3_5 <- ScheduledReceiptT_5[3]
r3_4 <- ScheduledReceiptT_4[3]
r3_3 <- ScheduledReceiptT_3[3]
r3_2 <- ScheduledReceiptT_2[3]
r3_1 <- ScheduledReceiptT_1[3]
r3_0 <- ScheduledReceipt[3]

r4_5 <- ScheduledReceiptT_5[4]
r4_4 <- ScheduledReceiptT_4[4]
r4_3 <- ScheduledReceiptT_3[4]
r4_2 <- ScheduledReceiptT_2[4]
r4_1 <- ScheduledReceiptT_1[4]
r4_0 <- ScheduledReceipt[4]

r5_5 <- ScheduledReceiptT_5[5]
r5_4 <- ScheduledReceiptT_4[5]
r5_3 <- ScheduledReceiptT_3[5]
r5_2 <- ScheduledReceiptT_2[5]
r5_1 <- ScheduledReceiptT_1[5]
r5_0 <- ScheduledReceipt[5]

r11 <- as.numeric()
r12 <- as.numeric()
r13 <- as.numeric()
r14 <- as.numeric()
r15 <- as.numeric()

r21 <- as.numeric()
r22 <- as.numeric()
r23 <- as.numeric()
r24 <- as.numeric()
r25 <- as.numeric()

r31 <- as.numeric()
r32 <- as.numeric()
r33 <- as.numeric()
r34 <- as.numeric()
r35 <- as.numeric()

r41 <- as.numeric()
r42 <- as.numeric()
r43 <- as.numeric()
r44 <- as.numeric()
r45 <- as.numeric()

r51 <- as.numeric()
r52 <- as.numeric()
r53 <- as.numeric()
r54 <- as.numeric()
r55 <- as.numeric()

# Inventory and Backlog
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)


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

#PERIOD 1
f.obj <- c(3,3,1,1,1,30,30)
f.con<- matrix(c(
  1,0,0,0,0,-1,0,
  0,1,0,0,0,0,-1,
  0,0,1,0,0,0,0,
  0,0,0,1,0,0,0,
  0,0,0,0,1,0,0,
                    
  0,0,0,0,0,1,0,
  0,0,0,0,0,0,1), nrow = 7, byrow = TRUE)

Constraint1.11 <- I10 - B10 - D10 + r1_1
Constraint1.21 <- I20 - B20 - D20 + r2_1
Constraint1.31 <- I30 - D30 + r3_3
Constraint1.41 <- I40 - D40 + r4_5
Constraint1.51 <- I50 - D50 + r5_1

Constraint2.11 <- B10 + D10
Constraint2.21 <- B20 + D20

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

f.dir <- c("=", "=", "=", "=", "=", "<=", "<=")

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

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]

#PERIOD 2
Constraint1.12 <- I11 - B11 - D11 + r1_0; Constraint1.12
Constraint1.22 <- I21 - B21 - D21 + r2_0; Constraint1.22
Constraint1.32 <- I31 - D31 + r3_2
Constraint1.42 <- I41 - D41 + r4_4  
Constraint1.52 <- I51 - D51 + r5_0

Constraint2.12 <- B11 + D11
Constraint2.22 <- B21 + D21

f2.rhs <- c(Constraint1.12, Constraint1.22, Constraint1.32, Constraint1.42,
            Constraint1.52, Constraint2.12, Constraint2.22)
f2.rhs

optimum2 <- lp(direction = "min", f.obj, f.con, f.dir, f2.rhs)
optimum2$solution

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]; B12
B22 <- optimum2$solution[7]


#PERIOD 3: works well when we have only B and I unknown (not r)
Constraint1.13 <- I12 - B12 - D12 + r11
Constraint1.23 <- I22 - B22 - D22 + r21
Constraint1.33 <- I32 - D32 + r3_1
Constraint1.43 <- I42 - D42 + r4_3 
Constraint1.53 <- I52 - D52 + r51

Constraint2.13 <- B12 + D12
Constraint2.23 <- B22 + D22

f3.rhs <- c(Constraint1.13, Constraint1.23, Constraint1.33, Constraint1.43,
            Constraint1.53, Constraint2.13, Constraint2.23)
f3.rhs
optimum3 <- lp(direction = "min", f.obj, f.con, f.dir, f3.rhs)
optimum3$solution

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]