main_illustrative_toy_problem.py

import argparse
from pathlib import Path
import logging

# import external packages
import matplotlib.pyplot as plt
import numpy as np
import pyomo.environ as pe
import pyomo.opt as po
logging.getLogger('pyomo.core').setLevel(logging.ERROR)

# import internal packages
from analysis import analysis
import dataio


output_path = Path("data/output/illustrative_toy_problem/")
output_path.mkdir(parents=True, exist_ok=True)

def get_nominal_data():
    return {
        'demand': (100, 100),
        'price_A': (1, 1),
        'price_B': (1.05, 1.05),
        'price_C': 2,
        'arc_inv_impact': 20,
        'arc_inv_cost': 1,
        'process_inv_impact': 20,
        'process_inv_cost': 2,
    }

def get_scenario_set(N=1000, random_seed=0):
    np.random.seed(random_seed)
    
    scenarios = np.random.uniform(size=(N,2))
    scenarios[:,0] += 0.5
    c1, c2 = scenarios[:,0], scenarios[:,1]
    scenarios[:,1] = c2*(1.5*c1 - 0.5*c1) + 0.5*c1

    return scenarios

def create_model(problem_instance, fixed_variable_info={}):    
    m = pe.ConcreteModel() 
    
    # initialize sets
    m.time_periods = pe.Set(initialize = [1,2])
    
    # initialize variables
    # continuous arc flow variables
    m.v_flow_SA_PA = pe.Var(m.time_periods, domain = pe.NonNegativeReals)
    m.v_flow_SB_PB = pe.Var(m.time_periods, domain = pe.NonNegativeReals)
    m.v_flow_PA_DC = pe.Var(m.time_periods, domain = pe.NonNegativeReals)
    m.v_flow_PB_DC = pe.Var(m.time_periods, domain = pe.NonNegativeReals)
    
    # integer investment variables
    m.v_inv_arc_SA_PA = pe.Var(m.time_periods, domain = pe.NonNegativeIntegers)
    m.v_inv_arc_SB_PB = pe.Var(m.time_periods, domain = pe.NonNegativeIntegers)
    m.v_inv_arc_PA_DC = pe.Var(m.time_periods, domain = pe.NonNegativeIntegers)
    m.v_inv_arc_PB_DC = pe.Var(m.time_periods, domain = pe.NonNegativeIntegers)
    m.v_inv_process_PA = pe.Var(m.time_periods, domain = pe.NonNegativeIntegers)
    m.v_inv_process_PB = pe.Var(m.time_periods, domain = pe.NonNegativeIntegers)

    # fix variables
    for var_name, fix_info in fixed_variable_info.items():
        var = getattr(m, var_name)
        for index, value in fix_info.items():
            var[index].fix(value)

    # initialize constraints
    @m.Constraint(m.time_periods)
    def satisfy_demand_C(m,t):
        return m.v_flow_PA_DC[t] + m.v_flow_PB_DC[t] == problem_instance['demand'][t-1]

    @m.Constraint(m.time_periods)
    def flow_balancing_PA(m,t):
        return m.v_flow_SA_PA[t] == m.v_flow_PA_DC[t]
    @m.Constraint(m.time_periods)
    def flow_balancing_PB(m,t):
        return m.v_flow_SB_PB[t] == m.v_flow_PB_DC[t]
        
    @m.Constraint(m.time_periods)
    def arc_capacity_SA_PA(m,t):
        return m.v_flow_SA_PA[t] <= sum(m.v_inv_arc_SA_PA[t2]
                                        for t2 in range(1,t+1)) * problem_instance['arc_inv_impact']
    @m.Constraint(m.time_periods)
    def arc_capacity_SB_PB(m,t):
        return m.v_flow_SB_PB[t] <= sum(m.v_inv_arc_SB_PB[t2]
                                        for t2 in range(1,t+1)) * problem_instance['arc_inv_impact']
    @m.Constraint(m.time_periods)
    def arc_capacity_PA_DC(m,t):
        return m.v_flow_PA_DC[t] <= sum(m.v_inv_arc_PA_DC[t2]
                                        for t2 in range(1,t+1)) * problem_instance['arc_inv_impact']
    @m.Constraint(m.time_periods)
    def arc_capacity_PB_DC(m,t):
        return m.v_flow_PB_DC[t] <= sum(m.v_inv_arc_PB_DC[t2]
                                        for t2 in range(1,t+1)) * problem_instance['arc_inv_impact']

    @m.Constraint(m.time_periods)
    def processing_capacity_PA(m,t):
        return m.v_flow_SA_PA[t] <= sum(m.v_inv_process_PA[t2]
                                        for t2 in range(1,t+1)) * problem_instance['process_inv_impact']
    @m.Constraint(m.time_periods)
    def processing_capacity_PB(m,t):
        return m.v_flow_SB_PB[t] <= sum(m.v_inv_process_PB[t2]
                                        for t2 in range(1,t+1)) * problem_instance['process_inv_impact']

    # initialize objective
    @m.Expression(m.time_periods)
    def revenue(m,t):
        return problem_instance['price_C'] * (m.v_flow_PA_DC[t] + m.v_flow_PB_DC[t])
    
    @m.Expression(m.time_periods)
    def supply_cost(m,t):
        return problem_instance['price_A'][t-1] * m.v_flow_SA_PA[t] + problem_instance['price_B'][t-1] * m.v_flow_SB_PB[t]
    
    @m.Expression(m.time_periods)
    def transport_investment_cost(m,t):
        return problem_instance['arc_inv_cost'] * (m.v_inv_arc_SA_PA[t] + m.v_inv_arc_SB_PB[t] + m.v_inv_arc_PA_DC[t] + m.v_inv_arc_PB_DC[t])

    @m.Expression(m.time_periods)
    def processing_investment_cost(m,t):
        return problem_instance['process_inv_cost'] * (m.v_inv_process_PA[t] + m.v_inv_process_PB[t])

    @m.Objective(sense=pe.maximize)
    def total_profit(m):
        return sum(m.revenue[t] - (m.supply_cost[t] + m.transport_investment_cost[t] + m.processing_investment_cost[t]) 
                   for t in m.time_periods)

    return m

def solve_model(m):
    solver = po.SolverFactory('gurobi')
    results = solver.solve(m)
    return results
    
def evaluate_feasibility(results):
    return results.solver.termination_condition != po.TerminationCondition.infeasible


def evaluate_objective(m):
    obj_value = pe.value(m.total_profit)
    return obj_value
    
def get_solution(m):
    solution = {
        str(var): {index: pe.value(var[index]) for index in var}
        for var in m.component_objects(pe.Var, active=True)
    }
    return solution

def alter_instance(problem_instance, scenario, decision_stages, current_stage):
    if decision_stages is None or current_stage == len(decision_stages)-1:
        problem_instance['price_A'] = (scenario[0], scenario[1])
    elif len(decision_stages) == 3 and current_stage == 1:
        problem_instance['price_A'] = (scenario[0], scenario[0])  # expect no deviation from t=1

    return problem_instance


def main(args):
    # First solve the nominal model:
    problem_instance = get_nominal_data()

    nominal_model = create_model(problem_instance)
    results = solve_model(nominal_model)
    feas_yn = evaluate_feasibility(results)
    if feas_yn:
        nominal_obj = evaluate_objective(nominal_model)
        nominal_solution = get_solution(nominal_model)

    random_seed = args.seed
    N = args.n
    scenarios = get_scenario_set(N=N, random_seed=random_seed)

    # sensitivity analysis results
    print("------------------------------------------------------------------------")
    print("SA")
    print("------------------------------------------------------------------------")
    SA_obj_path = output_path / f"tp_SA_obj_results_N={N}_new.csv"
    if SA_obj_path.exists() and not args.force_run:
        # Load results from .csv file
        SA_objective_results = dataio.read_array_results_from_csv(SA_obj_path)
        print("loaded prior SA results")
    else:
        decision_stages = None
        SA = analysis(decision_stages, scenarios, problem_instance,
                      create_model, solve_model, evaluate_feasibility,
                      evaluate_objective, get_solution, alter_instance,
                      analysis_method='SA', record_solutions=False,
                      verbose=args.verbose)

        SA_feasibility_results, SA_objective_results = SA.run(nominal_solution)
        # write results to .csv file
        dataio.write_array_results_to_csv(SA_obj_path, SA_objective_results)

    # robustness analysis in static setting
    print("------------------------------------------------------------------------")
    print("RA static")
    print("------------------------------------------------------------------------")
    RA_obj_path = output_path / f"tp_RA_obj_results_N={N}_new.csv"
    if RA_obj_path.exists() and not args.force_run:
        # Load results from .csv file
        RA_static_objective_results = dataio.read_array_results_from_csv(RA_obj_path)
        print("loaded prior RA statis results")
    else:
        decision_stages = {
            0: [('v_inv_arc_SA_PA',1),
                ('v_inv_arc_SB_PB',1),
                ('v_inv_arc_PA_DC',1),
                ('v_inv_arc_PB_DC',1),
                ('v_inv_process_PA',1),
                ('v_inv_process_PB',1),
                ('v_flow_SA_PA',1),
                ('v_flow_SB_PB',1),
                ('v_flow_PA_DC',1),
                ('v_flow_PB_DC',1),
                ('v_inv_arc_SA_PA',2),
                ('v_inv_arc_SB_PB',2),
                ('v_inv_arc_PA_DC',2),
                ('v_inv_arc_PB_DC',2),
                ('v_inv_process_PA',2),
                ('v_inv_process_PB',2),
                ('v_flow_SA_PA',2),
                ('v_flow_SB_PB',2),
                ('v_flow_PA_DC',2),
                ('v_flow_PB_DC',2)]
        }

        RA_static = analysis(decision_stages, scenarios, problem_instance,
                             create_model, solve_model, evaluate_feasibility,
                             evaluate_objective, get_solution, alter_instance,
                             analysis_method='RA', record_solutions=False,
                             verbose=args.verbose)

        RA_static_feasibility_results, RA_static_objective_results = RA_static.run(nominal_solution)

        # write results to .csv file
        dataio.write_array_results_to_csv(RA_obj_path, RA_static_objective_results)


    # robustness analysis in adaptive 3-stage setting
    print("------------------------------------------------------------------------")
    print("RA adaptive")
    print("------------------------------------------------------------------------")
    RA_adaptive_path = output_path / f"tp_RA_adaptive_obj_results_N={N}_new.csv"
    if RA_adaptive_path.exists() and not args.force_run:
        # Load results from .csv file
        RA_adaptive_objective_results = dataio.read_array_results_from_csv(RA_adaptive_path)
        print("loaded prior RA adaptive results")
    else:
        decision_stages = {
            0: [('v_inv_arc_SA_PA',1),
                ('v_inv_arc_SB_PB',1),
                ('v_inv_arc_PA_DC',1),
                ('v_inv_arc_PB_DC',1),
                ('v_inv_process_PA',1),
                ('v_inv_process_PB',1)],
            1: [('v_flow_SA_PA',1),
                ('v_flow_SB_PB',1),
                ('v_flow_PA_DC',1),
                ('v_flow_PB_DC',1),
                ('v_inv_arc_SA_PA',2),
                ('v_inv_arc_SB_PB',2),
                ('v_inv_arc_PA_DC',2),
                ('v_inv_arc_PB_DC',2),
                ('v_inv_process_PA',2),
                ('v_inv_process_PB',2)],
            2: [('v_flow_SA_PA',2),
                ('v_flow_SB_PB',2),
                ('v_flow_PA_DC',2),
                ('v_flow_PB_DC',2)]
        }

        RA_adaptive = analysis(decision_stages, scenarios, problem_instance,
                               create_model, solve_model, evaluate_feasibility,
                               evaluate_objective, get_solution, alter_instance,
                               analysis_method='RA', record_solutions=False,
                               verbose=args.verbose)

        RA_adaptive_feasibility_results, RA_adaptive_objective_results = RA_adaptive.run(nominal_solution)

        # write results to .csv file
        dataio.write_array_results_to_csv(RA_adaptive_path, RA_adaptive_objective_results)


    # plot results as histograms
    cmap = plt.get_cmap("tab10")

    dataio.plot_single_histogram("SA", args.plot_type, cmap(0), nominal_obj, SA_objective_results)
    dataio.plot_single_histogram("RA_static", args.plot_type, cmap(1), nominal_obj, RA_static_objective_results)
    dataio.plot_single_histogram("RA_adaptive", args.plot_type, cmap(2), nominal_obj, RA_adaptive_objective_results)


if __name__ == "__main__":
    parser = argparse.ArgumentParser()
    parser.add_argument("-f", "--force-run", action="store_true",
                        help="ignore any stored results and force (re)running (default: False)")
    parser.add_argument("-q", "--quiet", action="store_false", dest='verbose',
                        help="whether verbosity for analysis should be disabled (default False)")
    parser.add_argument("-n", type=int, default=1_000,
                        help="number of scenarios to run (default: 1000)")
    parser.add_argument("-s", "--seed", type=int, default=0,
                        help="random seed (default: 0)")
    parser.add_argument("--plot-type", type=str, default="pdf",
                        help="plot type, e.g. pdf, png, jpg (default: pdf)")
    args = parser.parse_args()
    main(args)