Notes

# -- coding: utf-8 --
"""
Spyder Editor

This is a temporary script file.
"""

# Importing necessary libraries
import pyomo.environ as pyo
from pyomo.opt import SolverFactory

# Input Data
components = ['Aluminum Frames', 'Carbon Fiber Frames', 'Manual Modules', 'Advanced Control Modules', 'Advanced Sensor Modules']
plants = ['Istanbul', 'Ankara', 'Izmir']
periods = [0, 1]

# Min and Max Demand for Components
minDemandComponents = [[0, 100, 200, 30, 100], [0, 100, 200, 30, 100], [0, 50, 100, 15, 100]]
maxDemandComponents = [[2000, 2000, 2000, 2000, 2000], [2000, 2000, 2000, 2000, 2000], [2000, 2000, 2000, 2000, 2000]]
minDemand = {}
maxDemand = {}
for i in range(len(components)):
for j in range(len(plants)):
minDemand[(components[i], plants[j])] = minDemandComponents[j][i]
maxDemand[(components[i], plants[j])] = maxDemandComponents[j][i]

# Min and Max Demand for Robotic Kits
minDemandRoboticKit = [0, 0, 0]
maxDemandRoboticKit = [200, 200, 200]
minDemandRK = {}
maxDemandRK = {}
for j in range(len(plants)):
minDemandRK[plants[j]] = minDemandRoboticKit[j]
maxDemandRK[plants[j]] = maxDemandRoboticKit[j]

# Production Cost for Components
productionCostComponents = [
[6, 19, 4, 10, 26],
[5, 18, 5, 11, 24],
[7, 20, 5, 12, 27]
]
productionCost = {}
for i in range(len(components)):
for j in range(len(plants)):
for k in range(len(periods)):
productionCost[(components[i], plants[j], periods[k])] = productionCostComponents[j][i] * (1.12 if k == 1 else 1)

# Production Cost for Robotic Kits
productionCostRoboticKits = [178, 175, 180]
productionCostRK = {}
for j in range(len(plants)):
for k in range(len(periods)):
productionCostRK[(plants[j], periods[k])] = productionCostRoboticKits[j] * (1.12 if k == 1 else 1)

# Inventory Cost for Components
inventoryCostComponents = {}
for i in range(len(components)):
for j in range(len(plants)):
for k in range(len(periods)):
inventoryCostComponents[(components[i], plants[j], periods[k])] = 0.08 * productionCost[(components[i], plants[j], periods[k])]

# Inventory Cost for Robotic Kits
inventoryCostRK = {}
for j in range(len(plants)):
for k in range(len(periods)):
inventoryCostRK[(plants[j], periods[k])] = 0.08 * productionCostRK[(plants[j], periods[k])]

# Selling Price for Components
sellingPriceComponents = [[10, 25, 8, 18, 40], [10, 25, 8, 18, 40], [12, 30, 10, 22, 45]]
sellingPrice = {}
for i in range(len(components)):
for j in range(len(plants)):
sellingPrice[(components[i], plants[j])] = sellingPriceComponents[j][i]

# Selling Price for Robotic Kits
sellingPriceRoboticKits = [290, 290, 310]
sellingPriceRK = {}
for j in range(len(plants)):
sellingPriceRK[plants[j]] = sellingPriceRoboticKits[j]

# Labor Data
laborData = [[1, 1.5, 1.5, 3, 4], [3.5, 3.5, 4.5, 4.5, 5], [3, 3.5, 4, 4.5, 5.5]]
labor = {}
for i in range(len(components)):
for j in range(len(plants)):
labor[(components[i], plants[j])] = laborData[j][i]

# Labor Availability
laborAvailabilityData = [12000, 15000, 22000]
laborAvailability = {}
for j in range(len(plants)):
laborAvailability[plants[j]] = laborAvailabilityData[j]

# Packing Data
packingData = [[4, 4, 5, 6, 6], [7, 7, 8, 9, 7], [7.5, 7.5, 8.5, 8.5, 8]]
packing = {}
for i in range(len(components)):
for j in range(len(plants)):
packing[(components[i], plants[j])] = packingData[j][i]

# Packing Availability
packingAvailabilityData = [20000, 40000, 35000]
packingAvailability = {}
for j in range(len(plants)):
packingAvailability[plants[j]] = packingAvailabilityData[j]

# Assembly Data
assemblyData = [65, 60, 65]
assembly = {}
for j in range(len(plants)):
assembly[plants[j]] = assemblyData[j]

# Assembly Availability
assemblyAvailabilityData = [5500, 5000, 6000]
assemblyAvailability = {}
for j in range(len(plants)):
assemblyAvailability[plants[j]] = assemblyAvailabilityData[j]

# Robotic Kit Requirements
req = [13, 13, 10, 3, 3]
roboticKitRequirements = {}
for i in range(len(components)):
roboticKitRequirements[components[i]] = req[i]

# Define the model
mdl = pyo.ConcreteModel()

# Sets
mdl.c = pyo.Set(initialize=components, doc='Components')
mdl.p = pyo.Set(initialize=plants, doc='Plants')
mdl.t = pyo.Set(initialize=periods, doc='Periods')

# Parameters
mdl.pMinDemand = pyo.Param(mdl.c, mdl.p, initialize=minDemand, doc='Minimum demand')
mdl.pMaxDemand = pyo.Param(mdl.c, mdl.p, initialize=maxDemand, doc='Maximum demand')
mdl.pMinDemandRK = pyo.Param(mdl.p, initialize=minDemandRK, doc='Minimum demand for robotic kits')
mdl.pMaxDemandRK = pyo.Param(mdl.p, initialize=maxDemandRK, doc='Maximum demand for robotic kits')
mdl.pLabor = pyo.Param(mdl.c, mdl.p, initialize=labor, doc='Labor requirement per unit')
mdl.pLaborAvailability = pyo.Param(mdl.p, initialize=laborAvailability, doc='Available labor minutes per plant')
mdl.pPacking = pyo.Param(mdl.c, mdl.p, initialize=packing, doc='Packing requirement per unit')
mdl.pPackingAvailability = pyo.Param(mdl.p, initialize=packingAvailability, doc='Available packing minutes per plant')
mdl.pAssembly = pyo.Param(mdl.p, initialize=assembly, doc='Assembly time per robotic kit')
mdl.pAssemblyAvailability = pyo.Param(mdl.p, initialize=assemblyAvailability, doc='Available assembly minutes per plant')
mdl.pRoboticKitRequirements = pyo.Param(mdl.c, initialize=roboticKitRequirements, doc='Component requirements per robotic kit')
mdl.pProductionCost = pyo.Param(mdl.c, mdl.p, mdl.t, initialize=productionCost, doc='Production cost per unit')
mdl.pProductionCostRK = pyo.Param(mdl.p, mdl.t, initialize=productionCostRK, doc='Production cost per robotic kit')
mdl.pSellingPrice = pyo.Param(mdl.c, mdl.p, initialize=sellingPrice, doc='Selling price per unit')
mdl.pSellingPriceRK = pyo.Param(mdl.p, initialize=sellingPriceRK, doc='Selling price per robotic kit')
mdl.pCarbonFiberLimit = pyo.Param(initialize=1000, doc='Total carbon fiber available per month')
mdl.pCarbonFiberPerFrame = pyo.Param(initialize=0.25, doc='Carbon fiber required per frame')
mdl.pInventoryCostRK = pyo.Param(mdl.p, mdl.t, initialize=inventoryCostRK, doc='Inventory cost for robotic kits')
mdl.pInventoryCostComponents = pyo.Param(mdl.c, mdl.p, mdl.t, initialize=inventoryCostComponents, doc='Inventory cost for components')


# Variables
mdl.vX_Produced = pyo.Var(mdl.c, mdl.p, mdl.t, bounds=(0.0, None), doc='Components produced')
mdl.vX_RK_Produced = pyo.Var(mdl.p, mdl.t, bounds=(0.0, None), doc='Robotic kits produced')
mdl.vInventory_Component = pyo.Var(mdl.c, mdl.p, mdl.t, bounds=(0.0, None), doc='Components sent to inventory')
mdl.vInventory_RK = pyo.Var(mdl.p, mdl.t, bounds=(0.0, None), doc='Robotic kits sent to inventory')
mdl.vSold_Component = pyo.Var(mdl.c, mdl.p, mdl.t, bounds=(0.0, None), doc='Components sold')
mdl.vSold_RK = pyo.Var(mdl.p, mdl.t, bounds=(0.0, None), doc='Robotic kits sold')

# Constraints

# 1. Labor Availability Constraint
def labor_availability_constraint(mdl, j, k):
return sum(mdl.vX_Produced[i, j, k] * mdl.pLabor[i, j] for i in mdl.c) - mdl.pLaborAvailability[j] <= 0
mdl.eLaborAvailability = pyo.Constraint(mdl.p, mdl.t, rule=labor_availability_constraint, doc="Labor availability constraint")

# 2. Packing Availability Constraint
def packing_availability_constraint(mdl, j, k):
return sum(mdl.vX_Produced[i, j, k] * mdl.pPacking[i, j] for i in mdl.c) - mdl.pPackingAvailability[j] <= 0
mdl.ePackingAvailability = pyo.Constraint(mdl.p, mdl.t, rule=packing_availability_constraint, doc="Packing availability constraint")

# 3. Carbon Fiber Constraint
def carbon_fiber_constraint(mdl, k):
return sum(mdl.vX_Produced['Carbon Fiber Frames', j, k] * mdl.pCarbonFiberPerFrame for j in mdl.p) - mdl.pCarbonFiberLimit <= 0
mdl.eCarbonFiber = pyo.Constraint(mdl.t, rule=carbon_fiber_constraint, doc="Carbon fiber usage constraint")

# 4. Assembly Availability Constraint
def assembly_availability_constraint(mdl, j, k):
return mdl.vX_RK_Produced[j, k] * mdl.pAssembly[j] - mdl.pAssemblyAvailability[j] <= 0
mdl.eAssemblyAvailability = pyo.Constraint(mdl.p, mdl.t, rule=assembly_availability_constraint, doc="Assembly availability constraint")

# 5. Minimum Demand Satisfaction (Components)
def min_demand_components_constraint(mdl, i, j, k):
return mdl.vSold_Component[i, j, k] - mdl.pMinDemand[i, j] >= 0
mdl.eMinDemandComponents = pyo.Constraint(mdl.c, mdl.p, mdl.t, rule=min_demand_components_constraint, doc="Minimum demand for components")

# 6. Maximum Demand Satisfaction (Sold Components)
def max_demand_components_constraint(mdl, i, j, k):
return mdl.vSold_Component[i, j, k] - mdl.pMaxDemand[i, j] <= 0
mdl.eMaxDemandComponents = pyo.Constraint(mdl.c, mdl.p, mdl.t, rule=max_demand_components_constraint, doc="Maximum demand for sold components")

# 7. Minimum Demand Satisfaction (Robotic Kits)
def min_demand_rk_constraint(mdl, j, k):
return mdl.vSold_RK[j, k] - mdl.pMinDemandRK[j] >= 0
mdl.eMinDemandRK = pyo.Constraint(mdl.p, mdl.t, rule=min_demand_rk_constraint, doc="Minimum demand for robotic kits")

# 8. Maximum Demand Satisfaction (Sold Robotic Kits)
def max_demand_rk_constraint(mdl, j, k):
return mdl.vSold_RK[j, k] - mdl.pMaxDemandRK[j] <= 0
mdl.eMaxDemandRK = pyo.Constraint(mdl.p, mdl.t, rule=max_demand_rk_constraint, doc="Maximum demand for robotic kits")

#9. Inventory Balance for Components (End of Month 1)
def inventory_balance_components_m1(mdl, i, j):
return mdl.vX_Produced[i, j, 0] - mdl.vSold_Component[i, j, 0] - mdl.vInventory_Component[i, j, 0] - mdl.vX_RK_Produced[j, 0]*mdl.pRoboticKitRequirements[i] == 0
mdl.eInventoryBalanceComponentsM1 = pyo.Constraint(mdl.c, mdl.p, rule=inventory_balance_components_m1, doc="Inventory balance for components at end of month 1")

# 10. Inventory Balance for Components (End of Month 2)
def inventory_balance_components_m2(mdl, i, j):
return mdl.vX_Produced[i, j, 1] + mdl.vInventory_Component[i, j, 0] - mdl.vSold_Component[i, j, 1] - mdl.vInventory_Component[i, j, 1] - mdl.vX_RK_Produced[j, 1]*mdl.pRoboticKitRequirements[i]== 0
mdl.eInventoryBalanceComponentsM2 = pyo.Constraint(mdl.c, mdl.p, rule=inventory_balance_components_m2, doc="Inventory balance for components at end of month 2")


# 11. Inventory Balance for Robotic Kits (End of Month 1)
def inventory_balance_rk_m1(mdl, j):
return mdl.vX_RK_Produced[j, 0] - mdl.vSold_RK[j, 0] - mdl.vInventory_RK[j, 0] == 0
mdl.eInventoryBalanceRKM1 = pyo.Constraint(mdl.p, rule=inventory_balance_rk_m1, doc="Inventory balance for robotic kits at end of month 1")

# 12. Inventory Balance for Robotic Kits (End of Month 2)
def inventory_balance_rk_m2(mdl, j):
return mdl.vX_RK_Produced[j, 1] + mdl.vInventory_RK[j, 0] - mdl.vSold_RK[j, 1] - mdl.vInventory_RK[j, 1] == 0
mdl.eInventoryBalanceRKM2 = pyo.Constraint(mdl.p, rule=inventory_balance_rk_m2, doc="Inventory balance for robotic kits at end of month 2")


# Objective Function
def objective_function(mdl):
# Sales revenue
revenue_components = sum(
mdl.vSold_Component[i, j, k] * mdl.pSellingPrice[i, j]
for i in mdl.c for j in mdl.p for k in mdl.t
)
revenue_rk = sum(
mdl.vSold_RK[j, k] * mdl.pSellingPriceRK[j]
for j in mdl.p for k in mdl.t
)

# Production costs
production_cost_components = sum(
mdl.vX_Produced[i, j, k] * mdl.pProductionCost[i, j, k]
for i in mdl.c for j in mdl.p for k in mdl.t
)
production_cost_rk = sum(
mdl.vX_RK_Produced[j, k] * mdl.pProductionCostRK[j, k]
for j in mdl.p for k in mdl.t
)

# Inventory costs
inventory_cost_components = sum(
mdl.vInventory_Component[i, j, k] * mdl.pInventoryCostComponents[i, j, k]
for i in mdl.c for j in mdl.p for k in mdl.t
)
inventory_cost_rk = sum(
mdl.vInventory_RK[j, k] * mdl.pInventoryCostRK[j, k]
for j in mdl.p for k in mdl.t
)

# Total profit
return (
revenue_components + revenue_rk
- production_cost_components - production_cost_rk
- inventory_cost_components - inventory_cost_rk
)

mdl.objective = pyo.Objective(rule=objective_function, sense=pyo.maximize, doc="Maximize profit")


# Solve the model
Solver = SolverFactory('glpk')
SolverResults = Solver.solve(mdl, tee=True)

# Display results
mdl.display()

# -- coding: utf-8 --
"""
Created on Wed Dec 6 12:43:14 2023

@author: Tolga
"""

# pyomo_sens_analysis.py

import numpy as np
import pandas as pd
import re

def read_LP_file(file_path_LP):
# Read the data into a DataFrame
df_LP = pd.read_csv(file_path_LP, header=0)

# Find the index of rows where the substrings in components_to_find exists in the df_LP's first column
components_to_find = ['c_l_', 'c_u_', 'c_e_', 'bound', 'end']
index_of_mdlLP_constraints=[]
index_of_mdlLP_variables=[]
for sub_comp in components_to_find:
if sub_comp != 'bound' and sub_comp != 'end':
index_of_mdlLP_constraints.append(df_LP[df_LP.iloc[:,0].str.contains(sub_comp)].index)
else:
index_of_mdlLP_variables.append(df_LP[df_LP.iloc[:,0].str.contains(sub_comp)].index)

LP_constraint_names = []
for index in index_of_mdlLP_constraints:
for sub_index in index:
LP_constraint_names.append(df_LP.iloc[sub_index,0][:-1])

df_LP_variables = df_LP.iloc[index_of_mdlLP_variables[0][0]+1:index_of_mdlLP_variables[1][0],0].copy()

# Split the entries in the 'Text' column based on the comma ('<=') delimiter
df_LP_variables_split = df_LP_variables.str.split('<=')

# Expand the list of split values into separate columns and store the second column in a list
LP_variable_names = df_LP_variables_split.apply(pd.Series).iloc[:,1].to_list()

return LP_variable_names, LP_constraint_names

def read_SA_file(file_path_SA):
# Specify the keyword to use for splitting
split_keyword = 'No.'

# Read the file and split data into chunks
with open(file_path_SA, 'r') as file:
data = file.read()

# Split the data into chunks using the specified keyword
chunks = data.split(split_keyword)

# Initialize a list to store DataFrames
data_frames = []

# Process each chunk
n=0
for chunk in chunks:
if not chunk.strip():
continue

# Split the chunk into lines
lines = chunk.strip().split('n')

# Identify the line containing column information
header_line = lines[2]

# Use regular expression to find positions of consecutive whitespaces
column_positions = [match.start() for match in re.finditer(r's+', header_line)]

# Initialize a list to store rows
rows = []

# Process each line in the chunk
for line in lines:
# Extract data based on identified column positions
data = [line[pos_start:pos_end].strip() for pos_start, pos_end in zip(column_positions, column_positions[1:]+[None])]
# Add the data to the current row
rows.append(data)

# Convert the list of rows into a DataFrame wrt being varible/constraint-related
if n>0:
if 'c_l_x' in rows[3][0] or 'c_u_x' in rows[3][0] or 'c_e_x' in rows[3][0] or 'c_e_ONE_VAR' in rows[3][0]:
df = pd.DataFrame(rows, columns=['Row name', 'St', 'Activity', ' Slack_Marginal', 'Lower bound_Upper bound',
'Activity range', 'Obj.Coeff range', 'Obj value at break point variable', 'Limiting variable'])
elif 'x' in rows[3][0] or 'ONE_VAR_CONSTANT' in rows[3][0]:
df = pd.DataFrame(rows, columns=['Column name', 'St', 'Activity', ' Obj coef_Marginal', 'Lower bound_Upper bound',
'Activity range', 'Obj.Coeff range', 'Obj value at break point variable', 'Limiting variable'])
else:
df = pd.DataFrame(rows)
n+=1


# Append the DataFrame to the list
data_frames.append(df)

#Concatenate dataframes wrt being varible/constraint-related
dfs_variables = pd.DataFrame(columns=['Column name', 'St', 'Activity', ' Obj coef_Marginal', 'Lower bound_Upper bound',
'Activity range', 'Obj.Coeff range', 'Obj value at break point variable', 'Limiting variable'])
dfs_constraints = pd.DataFrame(columns=['Row name', 'St', 'Activity', ' Slack_Marginal', 'Lower bound_Upper bound',
'Activity range', 'Obj.Coeff range', 'Obj value at break point variable', 'Limiting variable'])
n=0
for data_frame in data_frames:
if n>0:
if n<len(data_frames)-1:
df = data_frame[3:-5]
df = df.reset_index(drop=True)
elif n == len(data_frames)-1:
df = data_frame[3:-2]
df = df.reset_index(drop=True)
if 'Column name' in df.columns:
dfs_variables = pd.concat([dfs_variables, df])
elif 'Row name' in df.columns:
dfs_constraints = pd.concat([dfs_constraints, df])
n+=1

# Drop rows with empty strings in all columns
dfs_variables.replace('', pd.NA, inplace=True)
dfs_variables.dropna(how='all', inplace=True)
dfs_constraints.replace('', pd.NA, inplace=True)
dfs_constraints.dropna(how='all', inplace=True)

SA_constraint_names=dfs_constraints[dfs_constraints['Row name'].notna()]['Row name'].to_list()
SA_variable_names=dfs_variables[dfs_variables['Column name'].notna()]['Column name'].to_list()

return SA_variable_names, SA_constraint_names, dfs_variables, dfs_constraints


def reorganize_SA_report(file_path_SA, file_path_LP_labels, file_path_LP_nolabels):
"""
This is a brief explanation of the function "reorganize_SA_report" in the module "pyomo_sens_analysis".

Inputs:
- file_path_SA: Specify the path of your sensitivity analysis results file stored in a text file (e.g., 'D:YourFilePathYourFileName.txt').
- file_path_LP_labels: Specify the path of model file stored in a LP file with user-defined labels generated by "io_options={'symbolic_solver_labels': True}" (e.g., 'D:YourFilePathYourFileName.lp').
- file_path_LP_nolabels: Specify the path of model file stored in a LP file without user-defined labels generated by "io_options={'symbolic_solver_labels': False}" (e.g., 'D:YourFilePathYourFileName.lp').
Output:
Excel file containing sensitivity analysis results for obj. func. coeff. and constraint RHS in "Variables" and "Constraints" sheets, respectively.
"""
(SA_variable_names, SA_constraint_names, dfs_variables, dfs_constraints) = read_SA_file(file_path_SA)

(LP_variable_names_labels, LP_constraint_names_labels) =read_LP_file(file_path_LP_labels)

(LP_variable_names_nolabels, LP_constraint_names_nolabels) =read_LP_file(file_path_LP_nolabels)

LP_variable_name_mapping = {}
j=0
for i in LP_variable_names_nolabels:
LP_variable_name_mapping[i]=LP_variable_names_labels[j]
j+=1

LP_constraint_name_mapping = {}
j=0
for i in LP_constraint_names_nolabels:
LP_constraint_name_mapping[i]=LP_constraint_names_labels[j]
j+=1

dfs_variables_2 = dfs_variables.copy()
j=0
for i in np.arange(dfs_variables.shape[0]):
if j<len(LP_variable_name_mapping)-1 and isinstance(dfs_variables.iloc[i,0], str):
dfs_variables.iloc[i,0]=LP_variable_name_mapping[' ' + dfs_variables_2.iloc[i,0] + ' ']
j+=1

dfs_constraints_2 = dfs_constraints.copy()
j=0
for i in np.arange(dfs_constraints.shape[0]):
if j<len(LP_constraint_name_mapping)-1 and isinstance(dfs_constraints.iloc[i,0], str):
dfs_constraints.iloc[i,0]=LP_constraint_name_mapping[dfs_constraints_2.iloc[i,0]]
j+=1

# Create an ExcelWriter object
with pd.ExcelWriter('Sensitivity_Analysis_Report.xlsx') as writer:
# Write each DataFrame to a different sheet
dfs_variables.to_excel(writer, sheet_name='Variables',index=False)
dfs_constraints.to_excel(writer, sheet_name='Constraints',index=False)