Supplementary material to "Evaluating a hierarchical approach to landscape level harvest scheduling"

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

Section 1: List of variables which can be changed, modifying the size of the 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

Section 2: Adding packages


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

Section 3: Functions for the optimization models:

a) Monolithic model


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.

b) Department level optimizations:

Creates a set of 36 department level solutions which takes a small subset (6 solutions) of he NPI and PV trade-off for each of the 6 periods.


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

c) Integrated optimization problem:

Takes the department level solutions, and a selection of departments with complete information and solves. Through iterating multiple times, the solution will improve as the number of department level solutions increase.


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)

Section 4: Importing data


In [6]:
##Importing data:
os.chdir(path)
forest_data = pd.read_csv("data_heir_pseudo.dat",index_col=False, sep ="\t")

Section 5: Create data

The size of the problem depends on the number of departments, and the number of stands per department This problem follows the same logic as in the main test, 6 periods, first 4 are of 1 year long, last 2 are 3 years long.


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)

Section 6: Solving Optimization problems

a) Solving Monolithic problem


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))

b) Solving the Department level optimizations

This next section runs 36 optimizations for each department, striving to preliminarily cover the decision space. The focus is on both NPI and PV (according to the objective functions set earlier).


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))

c) Solving the integrated problem

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))