Fixed JuMP Code

###################Load Julia packages##########################################
using JuMP #package containing modeling features for optimization
using DataFrames #package facilitating input and output using dataframes
using CSV
using Cbc #load Cbc solver
using LinearAlgebra
using Statistics

###################Input data and initialize parameters##########################
forecastAndPrice=DataFrame!(CSV.File("RedTomatoForecastAndPrice_OffPeakDisc.csv"))
#forecastAndPrice=readtable("RedTomatoForecastAndPrice_OffPeakDisc.csv") #read forecast and price from file
demand = forecastAndPrice[:,2] #extract demand column of forecastAndPrice data frame
price  = forecastAndPrice[:,3] #extract price column of forecastAndPrice data frame

#Specify cost parameters
c_materials = 10 #Materials cost per unit
c_holding = 2 #Holding cost per unit/month
c_stockout = 5 #Stockout cost per unit/month
c_hiring = 300 #Hiring cost per worker
c_layoff = 500 #Layoff cost per worker
c_regular = 4 #Regular time wage per hour
c_overtime = 6 #Over time wage per hour
c_subcontract = 30 #Subcontracting per unit

#Specify other planning parameters
labReq = 4 #Labor hours per unit
hrsPerDay = 8 #hours per day
daysPerMonth = 20 #Days per month
maxOTPerMonth= 10 #Maximum over time per month
startingWorkForce= 80 #Maximum over time per month
endingWorkForce= 0 #ending workforce level
startingBackorder= 0 #starting backorder
endingBackorder= 0 #ending backorder
startInv= 1000 #starting inventory
endInv= 500 #ending inventory

########################Formulate optimization model####################################
nPeriods=size(demand,1) #obtain number of planning periods from forecast
#m = Model() #create Julia optimization model object
m = Model(with_optimizer(Cbc.Optimizer, logLevel=3))

#Create decision variables
@variable(m, W[0:nPeriods] >= 0,Int) #number of employees
@variable(m, H[1:nPeriods] >= 0, Int) #number of employees hired
@variable(m, L[1:nPeriods] >= 0, Int) #number of employess laid off
@variable(m, P[1:nPeriods] >= 0, Int) #number of units produced
@variable(m, Inv[0:nPeriods] >= 0) #number of units in inventory
@variable(m, S[0:nPeriods] >= 0) #number of units backordered
@variable(m, C[1:nPeriods] >= 0, Int) #number of units subcontracted
@variable(m, O[1:nPeriods] >= 0) #number of overtime hours

#Create objective function components
@expression(m, regTimeCost[t=1:nPeriods],  c_regular*hrsPerDay*daysPerMonth*W[t]) #regular time cost
@expression(m, OTCost[t=1:nPeriods],  c_overtime*O[t]) #over time cost
@expression(m, hiringCost[t=1:nPeriods],  c_hiring*H[t]) #hiring cost
@expression(m, layoffCost[t=1:nPeriods],  c_layoff*L[t]) #layoff cost
@expression(m, prodCost[t=1:nPeriods],  c_materials*P[t]) #production cost
@expression(m, invCost[t=1:nPeriods],  c_holding*Inv[t]) #inventory holding cost
@expression(m, stockOutCost[t=1:nPeriods],  c_stockout*S[t]) #stockout cost
@expression(m, subContCost[t=1:nPeriods], c_subcontract*C[t]) #subcontracting cost

#Add objective function components together and form "minimization" objective function
@objective(m, Min, sum(regTimeCost[t] + OTCost[t] + hiringCost[t]
+ layoffCost[t]  + prodCost[t] + invCost[t] + stockOutCost[t] + subContCost[t] for t=1:nPeriods))

#Create constraints
@constraints(m, begin
    const_initWorkforce, W[0]==startingWorkForce #set W[0] to initial workforce
    const_startInv, Inv[0] == startInv #set I[0] to initial inventory
    const_startBackorder,S[0]==startingBackorder #set S[0] to initial backorder

    const_workforceFlow[t=1:nPeriods], W[t] == W[t-1] + H[t] - L[t] #balance equations for workforce level
    const_invBal[t=1:nPeriods], Inv[t] == Inv[t-1] + P[t] + C[t] + S[t] - demand[t] - S[t-1] #balance equations for inventory

    const_prodCap[t=1:nPeriods], labReq*P[t] <= hrsPerDay*daysPerMonth*W[t] + O[t] #production time constraint
    const_overTime[t=1:nPeriods], O[t] <= maxOTPerMonth*W[t] #maximum overtime constraint

    const_endWorkforce, W[nPeriods] >= endingWorkForce #minimum ending workforce constraint
    const_endInv, Inv[nPeriods] >= endInv #minimum ending inventory constraint
    const_endBackorders, S[nPeriods] == endingBackorder #ending backorder constraint
end)

println("****Formulation****\n",m,"\n") #print optimization model

########################Solve model and output solution####################################
#status = solve(m) #solve optimization model
status = optimize!(m) #solve model

println("****Solve status and results****")
println("Optimization status: ", status,"\n") #print solver status
# println("Workforce: ", getvalue(W))

@show W
@show typeof(W)
@show getvalue.(W)
@show typeof(getvalue.(W))
W_=getindex.((getvalue.(W),), 1:nPeriods)
@show W_
@show typeof(W_)
H_=getvalue.(H)
L_=getvalue.(L)
P_=getvalue.(P)
C_=getvalue.(C)
S_=getindex.((getvalue.(S),), 1:nPeriods)
@show S_
@show typeof(S_)
Inv_=getindex.((getvalue.(Inv),), 1:nPeriods)
@show Inv_
@show typeof(Inv_)

#Create data frame with decision variables
decisionDF=DataFrame(W = W_, H = H_, L = L_, P = P_, C = C_, S = S_, Inv = Inv_)

#Create data frame with costs per period
costPerPeriodDF = DataFrame(regTimeCost=getvalue.(regTimeCost),OTCost=getvalue.(OTCost),
hiringCost=getvalue.(hiringCost),layoffCost=getvalue.(layoffCost),invCost=getvalue.(invCost),
stockOutCost=getvalue.(stockOutCost),prodCost=getvalue.(prodCost),subContCost=getvalue.(subContCost))

#Create data frame containing profit information
Revenue=dot(price,demand) #sum_{t = 1 to nperiods} price(t)*demand(t)
Cost=getobjectivevalue(m)
Profit=Revenue-Cost
profitDF=DataFrame(Revenue=Revenue,Cost=Cost,Profit=Profit)

#Process measures
Flowrate=mean(demand) #average demand satisfied by the aggregate plan
Inventory= 0.5*(getvalue.(Inv)[0] + getvalue.(Inv)[nPeriods]) + sum(getindex.((getvalue.(Inv),), 1:(nPeriods-1)))/nPeriods #average inventory
Flowtime=Inventory/Flowrate #flow time via little's law
processDF = DataFrame(Flowrate=Flowrate,Inventory=Inventory,Flowtime=Flowtime) #create data frame with process flow measures

println("****Optimal solution****\n",decisionDF,"\n") #print data frame containing decision variables created above
println(profitDF,"\n") #print profit
println("****Cost per period****\n",costPerPeriodDF,"\n") #print cost data frame created above
println("****Process flow measures****\n",processDF) #print data frame containing process flow measures

CSV.write("DecVar.csv", decisionDF) #write decision variable values to a csv file
CSV.write("profitSummary.csv", profitDF) #write profit information from data frame to a csv file
CSV.write("CostPerPeriod.csv", costPerPeriodDF) #write cost data frame to a csv file
CSV.write("ProcessFlowMeas.csv", processDF) #write process flow measures to a csv file