This iPython notebook contains code to reproduce the approach of solving the landscape level planning as detailed in the paper. The code is separated into several pieces, the sections of the code are listed at the end of this section. The problem is first solve as a monolithic problem (two times, once to calculate appropriate reference values for the second case), then solve a set of department level problems which span the trade-off between the net present income (NPI) and the present value (PV) of the stand at the end of the planning horizon.
The intent of this notebook is to provide a detailed example this approach, and highlight how changes in the size of the problem can impact the quality of the iterative approach and times taken to organize the data and solve the respective problems.
A note: The times taken to solve the entire integated model (including the set of department level solutions) in this iPython notebook are higher than the description listed in the manuscript. This is due to a lack of parallelization when solving the department level solutions.
Section 1: List of variables used throughout code
Section 2: Adding packages
Section 3: Optimization functions
a) Monolithic optimization function
b) Department level optimizations
c) Integerated hierarchical optimization function
Section 4: Importing data
Section 5: Creating structured data for optimization
Section 6: Solving optimization problems
a) Solving monolithic problem
b) Solving department level problems
c) Solving integrated problem
In [1]:
depts = 10 #The total number of departments used in the analysis
min_stands= 200 #The minimum number of stands per department
max_stands = 300 #The maximum number of stands per department
mono_time = 1800 #Number of seconds used to solve the monolithic problem
integrated_timelimit = 600 #Number of seconds used to solve each integrated problem
price_list = [50,30,50,30,40,20] #The expected price of the timber assortments [Spruce Logs, Spruce Pulp, Pine Logs, Pine Pulp, Birch Logs, Birch Pulp]
iters = 10 #The number of iterations used in the integrated problem
dept_used = 5 #The number of full department level information used in the integrated problem
prop_max_last_period = 0.6 #A proportion of the final period maximum yield
interest_rate = 0.03 #The interest rate used to calculate the NPI and PV.
Assortment_weights = [0.0133333]*36
Period_weights = [0.06,0.06,0.06,0.06,0.18,0.18]
path = "C:\MyTemp" #ADD PERSONAL PATH WHERE DATA IS LOCATED
min_area_stand = .1 #min area of stand in ha
max_area_stand = 6 #max area of stand in ha
In [2]:
#Add packages
import os
import pandas as pd
import random
from pyomo.environ import *
import itertools
from pyomo.opt import SolverFactory
import numpy as np
import pyomo
import datetime
In [3]:
# OPTIMIZATIONS
# FIRST: MONOLITHIC PROBLEM
# Creating decision variables, and initializign the problem
def monolithic(Whole_data,obj,reference_values):
times = [1,2,3,4,6,9]
mono_model = ConcreteModel()
mono_model.alternative= Var(set(Whole_data.index.values),bounds=(0,1),domain=NonNegativeReals)
mono_model.Periodic = Var(list(Whole_data)[3:-5]) #g_jt
mono_model.n_Periodic = Var(list(Whole_data)[3:-5],within=NonNegativeReals) #n_jt
mono_model.p_Periodic = Var(list(Whole_data)[3:-5],within=NonNegativeReals) #p_jt
mono_model.T_Periodic = Var(Periods) #g_t
mono_model.n_T_Periodic = Var(Periods,within=NonNegativeReals) #n_t
mono_model.p_T_Periodic = Var(Periods,within=NonNegativeReals) #p_t
mono_model.Dept_Decision = Var(set(range(0,len(Whole_data.Dept.unique())*len(Periods))),within=Binary) #first 1/6th P1, second 1/6th P2, ... last 1/6th P6
mono_model.NPI = Var()
mono_model.PV = Var()
mono_model.Objective = Var()
# Targets for the Goal programming portion of the problem. Initial values set high.
mono_model.reference = Param(range(len(list(Whole_data)[3:-5])),default = 10000000,mutable=True)
mono_model.component_reference = Param(Periods,default = 100000000,mutable=True)
if len(reference_values)>1:
for i in range(0,len(list(Whole_data)[3:-5])):
mono_model.reference[i] = reference_values[0][i]*prop_max_last_period
for j in range(1,7):
mono_model.component_reference[j] = reference_values[1][j-1]*prop_max_last_period
#CONSTRAINTS
#AREA CONSTRAINT
#For stand, a single schedule can be selected
for dept in Departments:
def onlyOneAlternativerule(model,stand):
return sum(model.alternative[idx]
for idx in Whole_data[(Whole_data["Dept"] == dept) & (Whole_data["Stand"] == stand)].index.values
)==1
setattr(mono_model,"alternative_const"+str(dept), Constraint(stands[dept],rule = onlyOneAlternativerule))
#Evaluating the periodic value flows of the timber assortments
def periodic_values(model,periodic_variable):
return model.Periodic[periodic_variable] == sum([model.alternative[idx]*Whole_data.at[idx,periodic_variable]*Whole_data.Area[idx]
for idx in Whole_data.index.values])
mono_model.periodic_values_const = Constraint(list(Whole_data)[3:-5],rule = periodic_values)
#Evaluating the periodic value flows of total timber flow
def T_periodic_values(model,T_variable):
return model.T_Periodic[T_variable] == sum([model.alternative[idx]*Whole_data.at[idx,periodic_variable]*Whole_data.Area[idx]
for idx in Whole_data.index.values for periodic_variable in list(Whole_data)[4+(T_variable-1)*7:10+(T_variable-1)*7]] )
mono_model.T_values_const = Constraint(Periods,rule = T_periodic_values)
#Periodic goal programming constraint for each timber assortment
def Periodic_GP(model,periodic_variable):
return model.reference[periodic_variable] == (model.Periodic[list(Whole_data)[periodic_variable+3]]-model.p_Periodic[list(Whole_data)[periodic_variable+3]]+model.n_Periodic[list(Whole_data)[periodic_variable+3]])
mono_model.Periodic_GP_constraint = Constraint(range(0,len(list(Whole_data)[3:-5])),rule=Periodic_GP)
#Periodic goal programming constraint for total timber flow
def Component_GP(model,T_variable):
return model.component_reference[T_variable] == (model.T_Periodic[T_variable]-model.p_T_Periodic[T_variable]+model.n_T_Periodic[T_variable])
mono_model.Component_GP_constraint = Constraint(Periods,rule=Component_GP)
#Present value accounting.
def mono_Present_value_rule(model):
return model.PV == sum([model.alternative[idx]*Whole_data.at[idx,PV_values]*Whole_data.Area[idx] for idx in Whole_data.index.values
for PV_values in list(Whole_data)[-5:-3]])/(1+interest_rate)**10
mono_model.Present_Value_rule_constraint = Constraint(rule=mono_Present_value_rule)
#Net Present income accounting.
def mono_Present_income_rule(model):
return model.NPI == sum([model.alternative[idx]*Whole_data.at[idx,NPI_Values]*Whole_data.Area[idx]*prices.at[0,NPI_Values]/(1+interest_rate)**time_t[NPI_Values]
for idx in Whole_data.index.values for NPI_Values in obj_list])
mono_model.Present_income_rule_constraint = Constraint(rule=mono_Present_income_rule)
#Constraining harvesting in the departments to only one time period
#for dept in Departments:
def Harv_dept(model,dept):
return sum(model.Dept_Decision[d2] for d2 in [dept+(i-1)*len(Departments) for i in Periods]) <= 1
mono_model.Harv_dept_constraint = Constraint(Departments,rule=Harv_dept)
#Calcuating H_dt
for period in Periods:
for dept in Departments:
def H_DTrule(model,stand):
return sum([model.alternative[idx] for idx in Whole_data[(Whole_data["Stand"] == stand) & (Whole_data["Dept"] == dept) & (Whole_data["Schedule"].isin([1+period,7+period,13+period]))].index.values])<= model.Dept_Decision[dept+(period-1)*len(Departments)] #added loop
setattr(mono_model,"H_DTrule_const"+str(dept)+"_"+str(period), Constraint(stands[dept],rule = H_DTrule))
#Setting the objective function
if obj == "APPROX_EVEN_FLOW":
mono_model.objective = Objective(expr=sum(
[mono_model.T_Periodic[T_variable]*T_periodic_weight[T_variable] for T_variable in Periods]+
[mono_model.Periodic[P_variable]*periodic_weight[P_variable] for P_variable in obj_list]+[mono_model.PV*0.0000001+mono_model.NPI*0.0000001]),sense=maximize)
elif obj == "GP":
mono_model.objective = Objective(expr=sum(
[mono_model.n_T_Periodic[T_variable]*T_periodic_weight[T_variable] +mono_model.p_T_Periodic[T_variable]*T_periodic_weight[T_variable] for T_variable in Periods]+
[mono_model.n_Periodic[P_variable]*periodic_weight[P_variable] + mono_model.p_Periodic[P_variable]*periodic_weight[P_variable] for P_variable in obj_list]
+[mono_model.PV*0.0000001+mono_model.NPI*0.0000001]),sense=minimize)
#Calling the optimization program
solver = "cplex"
opt = SolverFactory(solver)
opt.options['timelimit'] = mono_time # 30 minute time to solve.
opt.solve(mono_model,tee=True)
b1 = []
for a in obj_list:
b1.append(value(mono_model.Periodic[a]))
b1.append(value(mono_model.NPI))
b1.append(value(mono_model.PV))
print("Objective value: " + str(value(mono_model.objective)))
return b1, value(mono_model.objective) #returns the periodic flows of income, timber and the associated NPI and PV.
In [4]:
#This deletion of Hier_dataset could be removed once testing is done.
if "Hier_dataset" in locals():
del(Hier_dataset)
def pyomo_model(d,prices):
Periods = {1,2,3,4,5,6}
times = [1,2,3,4,6,9]
d1 = Department_data[Department_data["Dept"] == d]
d1 = d1.reset_index()
stands = d1["Stand"].unique().tolist()
dept_dataframe = pd.DataFrame(columns=obj_list+["NPI","PV"]) #creates a new dataframe that's empty
for t in Periods:
dept_model = ConcreteModel()
dept_model.alternative= Var(set(d1.index.values),bounds=(0,1),domain=NonNegativeReals)
dept_model.Periodic = Var(list(d1)[3:-5]) #g_jt
dept_model.Dept_Decision = Var(set(range(0,len(d1.Dept.unique())*len(Periods))),within=Binary) #first 1/6th P1, second 1/6th P2, ... last 1/6th P6
dept_model.NPI = Var()
dept_model.PV = Var()
dept_model.Objective = Var()
dept_model.lamb = Param(default = 0,mutable=True)
dept_model.PV_max = Param(default = 0,mutable=True)
dept_model.PV_min = Param(default = 0,mutable=True)
dept_model.t = Param(default = 1,mutable=True)
#CONSTRAINTS
#AREA CONSTRAINT
#For stand, a single schedule can be selected
def dept_onlyOneAlternativerule(model,stand):
return sum(model.alternative[idx]
for idx in d1[(d1["Stand"] == stand)].index.values
)==1
dept_model.alternative_const = Constraint(stands,rule = dept_onlyOneAlternativerule)
#Requires that all other periods except 't' do not have any harvests.
temp_periods = list(Periods)[:]
temp_periods.remove(t)
#Calcuating H_dt
def dept_only_t_harvest(model,periods):
return model.Dept_Decision[periods-1]== 0
dept_model.dept_only_t_harvest_const= Constraint(temp_periods,rule = dept_only_t_harvest)
#Evaluating the periodic value flows for all criteria of interest
def dept_periodic_values(model,periodic_variable):
return model.Periodic[periodic_variable] == sum([model.alternative[idx]*d1.at[idx,periodic_variable]*d1.Area[idx]
for idx in d1.index.values])
dept_model.periodic_values_const = Constraint(list(d1)[3:-5],rule = dept_periodic_values)
#Present value accounting.
def dept_Present_value_rule(model):
return model.PV == sum([model.alternative[idx]*d1.at[idx,PV_values]*d1.Area[idx] for idx in d1.index.values
for PV_values in list(d1)[-5:-3]])/(1+interest_rate)**10
dept_model.Present_Value_rule_constraint = Constraint(rule=dept_Present_value_rule)
#Net Present income accounting.
def dept_Present_income_rule(model):
return model.NPI == sum([model.alternative[idx]*d1.at[idx,NPI_Values]*d1.Area[idx]*prices.at[0,NPI_Values]/(1+interest_rate)**time_t[NPI_Values]
for idx in d1.index.values for NPI_Values in obj_list])
dept_model.Present_income_rule_constraint = Constraint(rule=dept_Present_income_rule)
#Constraining harvesting in the departments to only one time period
def dept_Harv_dept(model):
return sum(model.Dept_Decision[d2] for d2 in [(i-1) for i in Periods]) <= 1
dept_model.Harv_dept_constraint = Constraint(rule=dept_Harv_dept)
#Calcuating H_dt
for period in Periods:
def dept_H_DTrule(model,stand):
return sum([model.alternative[idx] for idx in d1[(d1["Stand"] == stand) & (d1["Dept"] == d) & (d1["Schedule"].isin([1+period,7+period,13+period]))].index.values])<= model.Dept_Decision[(period-1)] #added loop
setattr(dept_model,"H_DTrule_const"+"_"+str(period), Constraint(stands,rule = dept_H_DTrule))
#PV Constraint -- between Nadir and Ideal -- Nadir and Ideal calculated during first two optimizations.
def dept_PV_constraint_rule(model):
return model.PV >= model.PV_min + model.lamb*(model.PV_max -model.PV_min)
dept_model.PV_constraint_rule_constraint = Constraint(rule=dept_PV_constraint_rule)
#Parameters to change the problem, first is the objective function "NPI" maximized NPI,
#"PV" maximizes a combination of PV and NPI, where PV has a much higher weight
OBJ_RUN = ["NPI","PV","NPI","NPI","NPI","NPI"]
lamb = [0,0,0.2,0.4,0.6,0.8] #values of lambda
for k in range(0,len(OBJ_RUN)):
Obj = OBJ_RUN[k]
lamb_i = lamb[k]
if Obj == "NPI":
dept_model.objective = Objective(expr=dept_model.NPI,sense=maximize)
elif Obj == "PV":
dept_model.objective = Objective(expr=dept_model.PV*1000+dept_model.NPI/1000,sense=maximize)
dept_model.lamb=lamb_i
#RESTRICT TO SINGLE PERIOD
dept_model.preprocess()
solver = "cplex"
opt = SolverFactory(solver)
opt.solve(dept_model,tee=False)
#First two optimizations calculate the Nadir and Ideal values of PV
if k == 0:
dept_model.PV_min = dept_model.PV.value
elif k == 1:
dept_model.PV_max = dept_model.PV.value
dept_data = [value(dept_model.Periodic[i]) for i in obj_list]+[value(dept_model.NPI),value(dept_model.PV)]
ld = pd.DataFrame([dept_data], columns=obj_list+["NPI","PV"])
dept_dataframe = dept_dataframe.append(ld)
dept_model.del_component(dept_model.objective)
return dept_dataframe
In [5]:
def Integrated(D_DSET,H_DSET,ref_values,Departments_used,Departments_not_used):
j = 0
#Organize the full deparment level information for use in optimization
for i in Departments_not_used:
if j == 0:
DEPT_data = D_DSET[D_DSET["Dept"] == i]
j=j+1
elif j >= 1:
DEPT_data =DEPT_data.append(D_DSET[D_DSET["Dept"] == i])
DEPT_data = DEPT_data.reset_index()
sts = {}
for ii in Departments_not_used:
sts[ii] =DEPT_data[DEPT_data.Dept ==ii]["Stand"].unique()
H_DSET = H_DSET.reset_index()
H_DSET = H_DSET.drop('index', 1)
# Creating decision variables, and initializign the problem
Integrated_model = ConcreteModel()
Integrated_model.alternative= Var(set(H_DSET.index.values),bounds=(0,1),domain=Binary)
Integrated_model.FULL_alternative= Var(set(DEPT_data.index.values),bounds=(0,1),domain=NonNegativeReals)
Integrated_model.Periodic = Var(list(H_DSET)[0:-3]) #g_jt
Integrated_model.Periodic2 = Var(list(H_DSET)[0:-3]) #g_jt
Integrated_model.n_Periodic = Var(list(H_DSET)[0:-3],within=NonNegativeReals) #n_jt
Integrated_model.p_Periodic = Var(list(H_DSET)[0:-3],within=NonNegativeReals) #p_jt
Integrated_model.T_Periodic = Var(Periods) #g_t
Integrated_model.n_T_Periodic = Var(Periods,within=NonNegativeReals) #n_t
Integrated_model.p_T_Periodic = Var(Periods,within=NonNegativeReals) #p_t
Integrated_model.NPI = Var()
Integrated_model.PV = Var()
Integrated_model.Objective = Var()
Integrated_model.Dept_Decision = Var(set(range(0,len(DEPT_data.Dept.unique())*len(Periods))),within=Binary) #first 1/6th P1, second 1/6th P2, ... last 1/6th P6
# Targets for the Goal programming portion of the problem. Same as in the monolithic problem:
Integrated_model.reference = Param(range(len(list(H_DSET)[0:-3])),default = ref_values[0],mutable=True)
Integrated_model.component_reference = Param(Periods,default = 10000000,mutable=True)
for i in range(0,len(list(H_DSET)[0:-3])):
Integrated_model.reference[i] = ref_values[0][i]*prop_max_last_period
for j in Periods:
Integrated_model.component_reference[j] = ref_values[1][j-1]*prop_max_last_period
#CONSTRAINTS
#AREA CONSTRAINT
#For each department, one option can be selected
def Integrated_onlyOneAlternativeDeptrule(model,dept):
return sum(model.alternative[idx]
for idx in H_DSET[(H_DSET["Dept"] == dept)].index.values
)==1
Integrated_model.alternative_const_dept = Constraint(Departments_used,rule = Integrated_onlyOneAlternativeDeptrule)
#For those departments with full information, no option can be selected
def Integrated_onlyOneAlternativeDeptNOTUSEDrule(model,dept):
return sum(model.alternative[idx]
for idx in H_DSET[(H_DSET["Dept"] == dept)].index.values
)==0
Integrated_model.alternative_const_dept_not = Constraint(Departments_not_used,rule = Integrated_onlyOneAlternativeDeptNOTUSEDrule)
#FULL DEPARTMENT CONSTRAINTS
for dept in Departments_not_used:
def Integrated_onlyOneAlternativeSTANDrule(model,stand):
return sum(model.FULL_alternative[idx]
for idx in DEPT_data[(DEPT_data["Dept"] == dept) & (DEPT_data["Stand"] == stand)].index.values
)==1
setattr(Integrated_model,"alternative_STANDconst"+str(dept), Constraint(sts[dept],rule = Integrated_onlyOneAlternativeSTANDrule))
#Evaluating the periodic value flows for assorted timber values and income variables
def Integrated_periodic_values(model,periodic_variable):
return model.Periodic[periodic_variable] == sum(
[model.alternative[idx]*H_DSET.at[idx,periodic_variable] for idx in H_DSET.index.values]
+ [model.FULL_alternative[idx]*DEPT_data.at[idx,periodic_variable]*DEPT_data.Area[idx]
for idx in DEPT_data.index.values])
Integrated_model.periodic_values_const = Constraint(list(H_DSET)[0:-3],rule = Integrated_periodic_values)
#Periodic goal programming constraint for the assorted timber values and income variables.
def Integrated_Periodic_GP(model,periodic_variable):
return model.reference[periodic_variable] == (model.Periodic[list(H_DSET)[periodic_variable]]
-model.p_Periodic[list(H_DSET)[periodic_variable]]
+model.n_Periodic[list(H_DSET)[periodic_variable]])
Integrated_model.Periodic_GP_constraint = Constraint(range(0,len(list(H_DSET)[0:-3])),rule=Integrated_Periodic_GP)
#Evaluating the periodic value flows for the total timber flow
def Integrated_T_periodic_values(model,T_variable):
return model.T_Periodic[T_variable] == sum(
[model.alternative[idx]*H_DSET.at[idx,periodic_variable] for idx in H_DSET.index.values for periodic_variable in list(Hier_dataset)[(T_variable-1)*6:6+(T_variable-1)*6]]#+
+ [model.FULL_alternative[idx]*DEPT_data.at[idx,periodic_variable]*DEPT_data.Area[idx]
for idx in DEPT_data.index.values for periodic_variable in list(Hier_dataset)[(T_variable-1)*6:6+(T_variable-1)*6]]
)
Integrated_model.T_values_const = Constraint(Periods,rule = Integrated_T_periodic_values)
#Goal programming constraint for the total timber flows
def Integrated_Component_GP(model,T_variable):
return model.component_reference[T_variable] == (model.T_Periodic[T_variable]-model.p_T_Periodic[T_variable]+model.n_T_Periodic[T_variable])
Integrated_model.Component_GP_constraint = Constraint(Periods,rule=Integrated_Component_GP)
#Present value accounting. ******
def Integrated_Present_value_rule(model):
return model.PV == sum([model.FULL_alternative[idx]*DEPT_data.at[idx,PV_values]*DEPT_data.Area[idx] for idx in DEPT_data.index.values
for PV_values in list(DEPT_data)[-5:-3]]+
[model.alternative[idx]*H_DSET.at[idx,"PV"]for idx in H_DSET.index.values]
)/(1+interest_rate)**10
Integrated_model.Present_Value_rule_constraint = Constraint(rule=Integrated_Present_value_rule)
#Net Present income accounting.
def Integrated_Present_income_rule(model):
return model.NPI == sum([model.FULL_alternative[idx]*DEPT_data.at[idx,NPI_Values]*DEPT_data.Area[idx]*prices.at[0,NPI_Values]/(1+interest_rate)**time_t[NPI_Values]
for idx in DEPT_data.index.values for NPI_Values in obj_list]+
[model.alternative[idx]*H_DSET.at[idx,"NPI"]for idx in H_DSET.index.values])
Integrated_model.Present_income_rule_constraint = Constraint(rule=Integrated_Present_income_rule)
#Constraining harvesting in the departments to only one time period
#for dept in Departments:
def Integrated_Harv_dept(model,dept):
return sum(model.Dept_Decision[d2] for d2 in [dept+(i-1)*len(Departments_not_used) for i in Periods]) <= 1
Integrated_model.Harv_dept_constraint = Constraint(range(0,len(Departments_not_used)),rule=Integrated_Harv_dept)
#Calcuating H_dt
for period in Periods:
temp_value = -1
for dept in Departments_not_used:
temp_value = temp_value+1
def Integrated_H_DTrule(model,stand):
return sum([model.FULL_alternative[idx] for idx in DEPT_data[(DEPT_data["Stand"] == stand) &
(DEPT_data["Dept"] == dept) & (DEPT_data["Schedule"].isin([1+period,7+period,13+period]))].index.values])<= model.Dept_Decision[temp_value+(period-1)*len(Departments_not_used)] #added loop
setattr(Integrated_model,"H_DTrule_const"+str(dept)+"_"+str(period), Constraint(sts[dept],rule = Integrated_H_DTrule))
#The objective function
Integrated_model.objective = Objective(expr=sum([Integrated_model.n_T_Periodic[T_variable]*T_periodic_weight[T_variable] +Integrated_model.p_T_Periodic[T_variable]*T_periodic_weight[T_variable] for T_variable in Periods]+
[Integrated_model.n_Periodic[P_variable]*periodic_weight[P_variable] + Integrated_model.p_Periodic[P_variable]*periodic_weight[P_variable] for P_variable in obj_list]+
[Integrated_model.PV*0.0000001+Integrated_model.NPI*0.0000001]),sense=minimize)
solver = "cplex"
opt = SolverFactory(solver)
opt.options['timelimit'] = integrated_timelimit #A time limit in seconds to solve the problem
opt.solve(Integrated_model,tee=False)
#Create a list of periodic values provided by the solution
b1 = []
for a in obj_list:
b1.append(value(Integrated_model.Periodic[a]))
#Create department level specfic information on the periodic flows of the specific solution - to be integrated to the dataset
for DT in Departments_not_used:
a2 = []
b2 = []
for T_variable in range(1,7):
for tt in range(0,6):
periodic_variable = list(H_DSET)[(T_variable-1)*6:6+(T_variable-1)*6][tt]
a2.append(sum([value(Integrated_model.FULL_alternative[idx])*DEPT_data.at[idx,periodic_variable]*DEPT_data.Area[idx]
for idx in DEPT_data[DEPT_data.Dept == DT].index.values]))
ld = pd.DataFrame([a2+[0,0]+[DT]], columns=obj_list+["NPI","PV","Dept"])
H_DSET = H_DSET.append(ld)
return b1, H_DSET, value(Integrated_model.objective)
In [6]:
##Importing data:
os.chdir(path)
forest_data = pd.read_csv("data_heir_pseudo.dat",index_col=False, sep ="\t")
In [7]:
# Create data for hierarchical problem
#Create data for optimization, department data from the set of 1000 stands
for D in range(0,depts):
#random numbers used to select a variety of stands for the department
S = [random.randint(0,max(forest_data.Stand)) for x in range(random.randint(min_stands,max_stands))]
#Loop selecting the various stands from the dataset
for i in range(0,len(S)):
if i == 0:
fd = forest_data[forest_data.Stand == S[i]].assign(Area = random.uniform(min_area_stand,max_area_stand))
fd = fd.assign(Stand = i)
else:
s1 = forest_data[forest_data.Stand==S[i]].assign(Area = random.uniform(min_area_stand,max_area_stand))
s1 = s1.assign(Stand = i)
fd = fd.append(s1)
fd =fd.assign(Dept = D)
#Merging the department level data to a single dataset
if D == 0:
Department_data = fd
else:
Department_data = Department_data.append(fd)
Department_data = Department_data.reset_index()
#Create list of department names.
Departments = Department_data.Dept.unique().tolist()
#Create set of the number of periods:
Periods = {1,2,3,4,5,6}
#ORGANIZING DATA
stands = []
for i in Department_data.Dept.unique():
stands.append(set(Department_data[Department_data.Dept == i]["Stand"]) )
Departments = Department_data.Dept.unique().tolist()
dept_options = []
for i in Department_data.Dept.unique():
options = {}
for stand in stands[i]:
options[stand] =list(set(Department_data[(Department_data["Stand"] == stand)]["Schedule"]))
dept_options.append(options)
#Create list of objective values:
obj_list = []
for T_variable in Periods:
obj_list = obj_list + list(Department_data)[4+(T_variable-1)*7:10+(T_variable-1)*7]
#Used in NPI calculations
time_t = {}
time_values = [1]*6+[2]*6+[3]*6+[4]*6+[6]*6+[9]*6
for i in range(0,len(obj_list)):
time_t[str(obj_list[i])] = time_values[i]
#Creating a dictionary for the weights:
periodic_weight = {}
for i in range(0,len(obj_list)):
periodic_weight[str(obj_list[i])] = Assortment_weights[i]
T_periodic_weight = {}
for i in range(0,len(Periods)):
T_periodic_weight[i+1] = Period_weights[i]
#Setting the prices for each period to be the same.
prices = pd.DataFrame([price_list*6], columns=obj_list)
In [ ]:
#Use a proxy of the maximum sustainable flow as the targets for the goal programming portion of the problem.
#The model uses 80% of the reference value if GP model is used. As this case, the reference value is a theoretical max sustained yield
b1, mono_1_obj = monolithic(Department_data,"APPROX_EVEN_FLOW",[0])
reference_Values = []
reference_Values.append([0,b1[-8]/10,b1[-7]/10,b1[-6]/10,b1[-5]/10,b1[-4]/10,b1[-3]/10]*4+[0,3*b1[-8]/10,3*b1[-7]/10,3*b1[-6]/10,3*b1[-5]/10,3*b1[-4]/10,3*b1[-3]/10]*2)
reference_Values.append([sum(b1[-8:-2])/10]*4+[sum(b1[-8:-2])*3/10]*2)
tstart = datetime.datetime.now()
b2, mono_2_obj = monolithic(Department_data,"GP",reference_Values)
tend = datetime.datetime.now()
time =tend-tstart
print("Time taken to solve: " + str(time))
In [ ]:
#Starting time
tstart = datetime.datetime.now()
if "Hier_dataset" in locals():
del(Hier_dataset)
for dept in Departments:
print("Department: "+str(dept))
ld = pyomo_model(dept,prices)
ld.loc[:,'Dept'] = dept
if "Hier_dataset" in locals():
Hier_dataset = Hier_dataset.append(ld)
else:
Hier_dataset = ld
tend = datetime.datetime.now()
time =tend-tstart
print("Time in minutes: " +str(time.total_seconds()/60))
By changing the number of departments used, and the number of iterations run, the quality of the interated solution will change. Having too high number of departments used can impact the time taken to solve the problem. If there is a timelimit provided to the problem, an improved solution may not be found.
In [ ]:
#Re-organizing the reference values - the Income and PV values are not needed here.
R1 = [reference_Values[0][1:7]+reference_Values[0][1+7:7+7]+reference_Values[0][1+7*2:7+7*2]+reference_Values[0][1+7*3:7+7*3]
+reference_Values[0][1+7*4:7+7*4]+reference_Values[0][1+7*5:7+7*5]]+[[reference_Values[1][0]]+[reference_Values[1][1]]+
[reference_Values[1][2]]+[reference_Values[1][3]]+[reference_Values[1][4]]+[reference_Values[1][5]]]
#Loops through a pre-set number of iterations, key variables are the number of iterations (iters)
#and the number of full department information to be used (dept_used)
b = Hier_dataset
results= []
tstart = datetime.datetime.now()
for k in range(0,iters):
Departments_used = Departments[:]
Departments_not_used = random.sample(Departments_used,dept_used)
for i in Departments_not_used:
Departments_used.remove(i)
a1, b, c1 = Integrated(Department_data,b,R1,Departments_used,Departments_not_used)
results.append(c1)
print(str(k+1)+": "+str(c1))
tend = datetime.datetime.now()
time =tend-tstart
print()
print("Time taken to solve: " + str(time))
print()
print("Objective function of each iteration: ")
for r in range(0,len(results)):
print(str(r+1)+": "+str(results[r]))
print()
print("Objective function of the monolithic problem: " + str(mono_2_obj))