HW2Problem3Code.py

# -*- coding: utf-8 -*-
"""
Created on Sat Mar 20 09:50:25 2021

@author: Moon
"""

# Import of the pyomo module
from pyomo.environ import *
import numpy as np
import math
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

# Creation of a Concrete Model
Optimal_values_I=[]
result_w=[]
Optimal_values=[]
M=100000

B=75000

model_I = ConcreteModel()

# Initialize state and actions
model_I.i = Set(initialize=['1','2','3'], doc='State i')
model_I.j = Set(initialize=['1','2','3'], doc='State j')
model_I.a = Set(initialize=['a1','a2', 'a3'], doc='Action')

# User cost
model_I.u_cost = Param(model_I.i, initialize={'1':10,'2':20,'3':30}, doc='UserCost')

# Action costs
d_cost = {
    ('a1','1'): 0.0,
    ('a1','2'): 0.0,
    ('a1','3'): 0.0,
    ('a2','1'): 5.0,
    ('a2','2'): 10.0,
    ('a2','3'): 15.0,
    ('a3','1'): 20.0,
    ('a3','2'): 20.0,
    ('a3','3'): 20.0,
    }

model_I.a_cost = Param(model_I.a, model_I.i, initialize=d_cost, doc='ActionCost')

# Transition Propability

P_DN={
    ('1','1') : 0.5,
    ('1','2') : 0.5,
    ('1','3') : 0,
    ('2','1') : 0,
    ('2','2') : 0.5,
    ('2','3') : 0.5,
    ('3','1') : 0,
    ('3','2') : 0,
    ('3','3') : 1,
    }

P_RM={
    ('1','1'):0.8,
    ('1','2'):0.2,
    ('1','3'):0,
    ('2','1'):0,
    ('2','2'):0.8,
    ('2','3'):0.2,
    ('3','1'):0,
    ('3','2'):0,
    ('3','3'):1,
    }
P_Resurf={
    ('1','1'): 1.0,
    ('1','2'): 0,
    ('1','3'): 0,
    ('2','1'): 0.8,
    ('2','2'): 0.2,
    ('2','3'): 0,
    ('3','1'): 0.6,
    ('3','2'): 0.4,
    ('3','3'): 0,
    }

P={"a1":P_DN,"a2":P_RM,"a3":P_Resurf}

model_I.P = Param(model_I.a, initialize=P, doc='Transition probability')
model_I.w = Var(model_I.a, model_I.i, bounds=(0,1), doc='% of the pavement with certain action')

# Define constraints
def const_1(model_I,a,i):
    return model_I.w[a,i] >= 0

model_I.const_1 = Constraint(model_I.a, model_I.i, rule=const_1, doc='Wai>=0')
   
def const_2(model_I):
    sum_wai2 = 0.0
    for ka in model_I.a:
        for ki in model_I.i:
            sum_wai2 += model_I.w[ka,ki]
    
    return sum_wai2==1

model_I.const_2 = Constraint(rule=const_2,doc='w_constraints = 1')
model_I.const_2.display()

def const_3(model_I,i):
    sum_wai3=0
    sum_pw3=0
    for ka in model_I.a:
        sum_wai3+=model_I.w[ka,i]
        for kj in model_I.j:
            sum_pw3 += model_I.P[ka][kj,i]*model_I.w[ka,kj]
 
    return sum_wai3 == sum_pw3

model_I.const_3 = Constraint(model_I.i, rule=const_3,doc='sum_wai==sum_pw')
   
def const_4(model_I):
    sum_cw = 0.0
    for ka in model_I.a:
        for ki in model_I.i:
            sum_cw += model_I.w[ka,ki]*model_I.a_cost[ka,ki]
    
    return sum_cw*M<=B

model_I.const_4 = Constraint(rule=const_4,doc='budget constraints ')
def const_5(model_I,i):
    sum_wai=0
    for ka in model_I.a:
        sum_wai+=model_I.w[ka,i]
    
    return sum_wai >= 0.0
model_I.const_5 = Constraint(model_I.i,rule=const_5,doc='w_constraints >= 0')

def const_6(model_I,i):
    sum_wai=0
    for ka in model_I.a:
        sum_wai+=model_I.w[ka,i]
    return sum_wai <= 1

model_I.const_6 = Constraint(model_I.i,rule=const_6,doc='w_constraints <= 1.0')
def objective_rule(model_I):
    obj=0
    for ki in model_I.i:
        for ka in model_I.a:
            obj+=model_I.w[ka,ki]*(model_I.u_cost[ki] + model_I.a_cost[ka,ki])

    return obj*M

model_I.objective = Objective(rule=objective_rule, sense=minimize, doc='Define objective function')


def pyomo_postprocess(options=None, instance=None, results=None):
  model_I.w.display()

# This is an optional code path that allows the script to be run outside of
# pyomo command-line.  For example:  python transport.py
if __name__ == '__main__':
    import pyomo.environ
    from pyomo.opt import SolverFactory
    opt = SolverFactory("glpk")
    results = opt.solve(model_I)
    
    # Printout results
    results.write()
    print("\nDisplaying Solution\n" + '-'*60)
    pyomo_postprocess(None, model_I, results)

print(model_I.objective())
  

# Creation of a Concrete Model
for T in range (25,226,25):
    model_F = ConcreteModel()
    
    T0=T-1
    model_F.i = Set(initialize=['1','2','3'], doc='i')
    model_F.j = Set(initialize=['1','2','3'], doc='j')
    model_F.a = Set(initialize=['a1','a2', 'a3'], doc='a')
    model_F.t = Set(initialize=np.arange(T)+1, doc='states')
    model_F.t0 = Set(initialize=np.arange(T0)+1, doc='states')
    model_F.u_cost = Param(model_F.i, initialize={'1':10,'2':20,'3':30}, doc='user_cost')
    model_F.a_cost = Param(model_F.a, model_F.i, initialize=d_cost, doc='action_cost[s,a]')
    model_F.P = Param(model_F.a, initialize=P, doc='tr_probability_routine_maintenance P_ji')
    model_F.w = Var(model_F.a, model_F.i,model_F.t, bounds=(0,1), doc='% of the pavement with certain action')

    def const_1(model,a,i,t):
        #print(model.w[a,i,t] >= 0)
        return model.w[a,i,t] >= 0


    model_F.const_1 = Constraint(model_F.a, model_F.i, model_F.t, rule=const_1, doc='Wai>=0')


    def const_2(model_F,t):
        sum_wai2 = 0.0
        for ka in model_F.a:
            for ki in model_F.i:
                sum_wai2 += model_F.w[ka,ki,t]
        #print(sum_wai2==1)
        return sum_wai2==1.0
    model_F.const_2 = Constraint(model_F.t,rule=const_2,doc='sum_wai2==1')


    def const_3(model_F,i,t0):
        sum_wai3=0
        sum_pw3=0
        #t0=[kt+1 for kt in t1]
        for ka in model_F.a:
            sum_wai3+=model_F.w[ka,i,t0+1]
            for kj in model_F.j:
                sum_pw3 += model_F.P[ka][kj,i]*model_F.w[ka,kj,t0]
        return sum_wai3 == sum_pw3

    model_F.const_3 = Constraint(model_F.i,model_F.t0,rule=const_3,doc='sum_wai==sum_pw')

    def const_4(model_F,t):
        sum_cw = 0.0
        for ka in model_F.a:
            for ki in model_F.i:
                sum_cw += model_F.w[ka,ki,t]*model_F.a_cost[ka,ki]
        #print(sum_cw*M <=B)
        return sum_cw*M <= B

    model_F.const_4 = Constraint(model_F.t,rule=const_4,doc='budget constraints')

    def const_5(model_F,i,t):
        sum_wai5=0
        for ka in model_F.a:
            sum_wai5+=model_F.w[ka,i,t]
        #print(sum_wai5>= 0)
        return sum_wai5 >= 0.0

    model_F.const_5 = Constraint(model_F.i,model_F.t,rule=const_5,doc='w_constraints >= 0')
    #model.const_5.display()

    def const_6(model_F,i,t):
        sum_wai6=0
        for ka in model_F.a:
            sum_wai6+=model_F.w[ka,i,t]
        #print(sum_wai6<= 1)
        return sum_wai6 <= 1

    model_F.const_6 = Constraint(model_F.i,model_F.t,rule=const_6,doc='w_constraints <= 1.0')

    def const_7(model_F,i):
        sum_wai7=0
        for ka in model_F.a:
            sum_wai7+=model_F.w[ka,i,1]
        if i == '1':
            return sum_wai7 == 0
        elif i == '2':
            return sum_wai7 == 1
        elif i == '3':
            return sum_wai7 == 0

    model_F.const_7 = Constraint(model_F.i,rule=const_7,doc='boundary1')


    w_t={
    ('a1', '1'):0.0,
    ('a1', '2'):0.043,
    ('a1', '3'):0.872,
    ('a2', '1'):0.064,
    ('a2', '2'):0.0,
    ('a2', '3'):0.0,
    ('a3', '1'):0.0,
    ('a3', '2'):0.0,
    ('a3', '3'):0.021,
    }

    def const_8(model_F,i):
        sum_wai8=0
        sum_pw8=0
        for ka in model_F.a:
            sum_wai8+=model_F.w[ka,i,T]
            sum_pw8 += value(model_I.w[ka,i])
        #print(sum_wai8 ==sum_pw8)
        return sum_wai8 == sum_pw8

    model_F.const_8 = Constraint(model_F.i, rule=const_8, doc='boundary2')


    def objective_rule(model_F):
        obj=0
        for kt in model_F.t:
            for ka in model_F.a:
                for ki in model_F.i:
                    obj+=(model_F.w[ka,ki,kt])*(model_F.u_cost[ki] + model_F.a_cost[ka,ki])*(0.97)**(kt)
        return obj*M + (0.97**T/(1-0.97))*model_I.objective()

    model_F.objective = Objective(rule=objective_rule, sense=minimize, doc='Define objective function')

    ## Display of the output ##
    # Display x.l, x.m ;
    def pyomo_postprocess(options=None, instance=None, results=None):
      model_F.w.display()


    try:

        if __name__ == '__main__':
            # This emulates what the pyomo command-line tools does
            import pyomo.environ
            from pyomo.opt import SolverFactory
            opt = SolverFactory("glpk",validate=False)
            results = opt.solve(model_F)
            #sends results to stdout
            results.write()
            print("\nDisplaying Solution\n" + '-'*60)
            pyomo_postprocess(None, model_F, results)

            objective_value_F = value(model_F.objective)
            objective_value_I = value(model_I.objective)
    except:
            objective_value_F= np.NAN
            objective_value_I= np.NAN


    Optimal_values.append([T, B, objective_value_F, objective_value_I])



df_optimal=pd.DataFrame(Optimal_values,columns=["T","B","ObjectiveValueF", "ObjectiveValueI"])

df_optimal.to_csv("df_optimal.csv")


# Plot for different budget and optimal costs

import matplotlib.pyplot as plt


    
fig,ax=plt.subplots(figsize=(10,6))

df_plot=df_optimal[df_optimal['B']==B]

plt.plot(df_plot["T"],df_plot['ObjectiveValueF'])  

plt.xlabel('Short Term Horizon in Years', size="16")
plt.ylabel('Optimal Cost',size="16")
plt.title('Finite Horizon Problem, Budget :'+str(B),size="16" )
plt.tight_layout()

plt.savefig(str(B)+"variation.png",dpi=500)
plt.show()