|
|
import pandas as pd |
|
|
import numpy as np |
|
|
import random |
|
|
import copy |
|
|
import optuna |
|
|
import matplotlib.pyplot as plt |
|
|
import seaborn as sns |
|
|
from skopt import gp_minimize |
|
|
from skopt.space import Real, Categorical |
|
|
from skopt.utils import use_named_args |
|
|
from statsmodels.tsa.arima.model import ARIMA |
|
|
import numpy_financial as npf |
|
|
|
|
|
optuna.logging.set_verbosity(optuna.logging.WARNING) |
|
|
|
|
|
PROJECT_YEARS = 15 |
|
|
BASE_CAPACITY_KTA = 300 |
|
|
INFLATION_RATE = 0.015 |
|
|
TAX_RATE = 0.10 |
|
|
DEPRECIATION_YEARS = 15 |
|
|
|
|
|
TECHNOLOGY_DATA = { |
|
|
"JNC": {"capex_base_M": 180.0, "opex_base_cents_kg": 150.0}, |
|
|
"Hoechst_AG": {"capex_base_M": 230.0, "opex_base_cents_kg": 155.0}, |
|
|
"BF_Goodrich": {"capex_base_M": 280.0, "opex_base_cents_kg": 160.0}, |
|
|
"Shin_Etsu_1991": {"capex_base_M": 260.0, "opex_base_cents_kg": 155.0}, |
|
|
"Shin_Etsu_2004": {"capex_base_M": 1500.0, "opex_base_cents_kg": 150.0}, |
|
|
"Vinnolit": {"capex_base_M": 240.0, "opex_base_cents_kg": 155.0}, |
|
|
"QVC_Qatar": {"capex_base_M": 200.0, "opex_base_cents_kg": 145.0}, |
|
|
"SP_Chemicals": {"capex_base_M": 250.0, "opex_base_cents_kg": 145.0}, |
|
|
"Engro_Pakistan": {"capex_base_M": 1400.0, "opex_base_cents_kg": 155.0}, |
|
|
"Formosa_BR_USA": {"capex_base_M": 380.0, "opex_base_cents_kg": 150.0}, |
|
|
"Shintech_USA_Exp": {"capex_base_M": 1700.0, "opex_base_cents_kg": 160.0}, |
|
|
"Zhongtai_China": {"capex_base_M": 1100.0, "opex_base_cents_kg": 155.0}, |
|
|
"Shintech_USA_2021": {"capex_base_M": 1700.0, "opex_base_cents_kg": 160.0}, |
|
|
"Reliance_India_2024": {"capex_base_M": 2200.0, "opex_base_cents_kg": 155.0}, |
|
|
"Orbia_Germany_2023": {"capex_base_M": 180.0, "opex_base_cents_kg": 155.0}, |
|
|
"Westlake_USA_2022": {"capex_base_M": 850.0, "opex_base_cents_kg": 160.0} |
|
|
} |
|
|
|
|
|
STRATEGY_DATA = { |
|
|
'Integrated_Production': {'sourcing_cost_per_ton_pvc': 450.0, 'byproducts': {'caustic_soda_ton': 1.1, 'surplus_edc_ton': 0.523}}, |
|
|
'Purchase_VCM': {'sourcing_cost_per_ton_pvc': 650.0, 'byproducts': {'caustic_soda_ton': 0, 'surplus_edc_ton': 0}} |
|
|
} |
|
|
|
|
|
PRODUCT_PRICES_USD_PER_TON = { |
|
|
'pvc_s65_export': 1100, 'pvc_s70_domestic': 950, |
|
|
'byproduct_caustic_soda': 450, 'byproduct_surplus_edc': 170 |
|
|
} |
|
|
|
|
|
OPTIMIZATION_SPACE = { |
|
|
'capacity_kta': (500, 600), |
|
|
'technology': ["Engro_Pakistan", "Shin_Etsu_2004"], |
|
|
'sourcing_strategy': ['Integrated_Production'], |
|
|
'export_market_mix': (0.6, 0.8), |
|
|
'sell_byproducts': [True] |
|
|
} |
|
|
|
|
|
def forecast_prices(base_price, years=PROJECT_YEARS, cagr=0.04): |
|
|
historical = [base_price * (1 + cagr + random.uniform(-0.03, 0.03)) ** i for i in range(-5, 0)] |
|
|
model = ARIMA(historical, order=(1, 1, 0)) |
|
|
fit = model.fit() |
|
|
forecast = fit.forecast(steps=years) |
|
|
return [p * (1 + INFLATION_RATE) ** i for i, p in enumerate(forecast)] |
|
|
|
|
|
def calculate_project_kpis(**params): |
|
|
tech_data = TECHNOLOGY_DATA[params['technology']] |
|
|
scaling_factor = (params['capacity_kta'] / BASE_CAPACITY_KTA) ** 0.65 |
|
|
total_capex = tech_data['capex_base_M'] * scaling_factor * 1_000_000 |
|
|
capacity_tons = params['capacity_kta'] * 1000 |
|
|
|
|
|
price_s65_forecast = forecast_prices(PRODUCT_PRICES_USD_PER_TON['pvc_s65_export']) |
|
|
price_s70_forecast = forecast_prices(PRODUCT_PRICES_USD_PER_TON['pvc_s70_domestic']) |
|
|
byproduct_caustic_forecast = forecast_prices(PRODUCT_PRICES_USD_PER_TON['byproduct_caustic_soda']) |
|
|
byproduct_edc_forecast = forecast_prices(PRODUCT_PRICES_USD_PER_TON['byproduct_surplus_edc']) |
|
|
|
|
|
cash_flows = [-total_capex] |
|
|
depreciation_annual = total_capex / DEPRECIATION_YEARS |
|
|
|
|
|
for year in range(1, PROJECT_YEARS + 1): |
|
|
infl_factor = (1 + INFLATION_RATE) ** (year - 1) |
|
|
price_s65 = price_s65_forecast[year - 1] |
|
|
price_s70 = price_s70_forecast[year - 1] * 0.95 if params['capacity_kta'] >= 550 else price_s70_forecast[year - 1] |
|
|
revenue_export = (capacity_tons * params['export_market_mix']) * price_s65 |
|
|
revenue_domestic = (capacity_tons * (1 - params['export_market_mix'])) * price_s70 |
|
|
pvc_revenue = revenue_export + revenue_domestic |
|
|
|
|
|
byproduct_revenue = 0 |
|
|
if params['sourcing_strategy'] == 'Integrated_Production' and params['sell_byproducts']: |
|
|
byproducts = STRATEGY_DATA['Integrated_Production']['byproducts'] |
|
|
byproduct_revenue += (byproducts['caustic_soda_ton'] * capacity_tons * byproduct_caustic_forecast[year - 1]) |
|
|
byproduct_revenue += (byproducts['surplus_edc_ton'] * capacity_tons * byproduct_edc_forecast[year - 1]) |
|
|
|
|
|
total_revenue = pvc_revenue + byproduct_revenue |
|
|
|
|
|
opex_sourcing = STRATEGY_DATA['Integrated_Production']['sourcing_cost_per_ton_pvc'] * capacity_tons * infl_factor |
|
|
opex_base = (tech_data['opex_base_cents_kg'] / 100) * capacity_tons * infl_factor |
|
|
total_opex = opex_sourcing + opex_base |
|
|
|
|
|
ebitda = total_revenue - total_opex |
|
|
taxable_income = ebitda - depreciation_annual |
|
|
taxes = max(taxable_income * TAX_RATE, 0) |
|
|
net_income = taxable_income - taxes |
|
|
free_cash_flow = net_income + depreciation_annual |
|
|
cash_flows.append(free_cash_flow) |
|
|
|
|
|
irr = npf.irr(cash_flows) * 100 if npf.irr(cash_flows) > 0 else -1.0 |
|
|
cumulative_cash_flow = np.cumsum(cash_flows) |
|
|
payback_period_years = np.where(cumulative_cash_flow > 0)[0] |
|
|
payback_period = payback_period_years[0] + 1 if len(payback_period_years) > 0 else float('inf') |
|
|
annual_profit = np.mean([cf for cf in cash_flows[1:] if cf > 0]) |
|
|
|
|
|
return {"irr": irr, "annual_profit": annual_profit, "total_capex": total_capex, "payback_period": payback_period} |
|
|
|
|
|
def run_optimizations_without_ml(): |
|
|
print("\n--- Running Optimization Algorithms ---") |
|
|
results = [] |
|
|
results.append(run_bayesian_optimization()) |
|
|
results.append(run_genetic_algorithm()) |
|
|
results.append(run_optuna_direct()) |
|
|
return results |
|
|
|
|
|
def run_genetic_algorithm(): |
|
|
print("Running Genetic Algorithm (GA)...") |
|
|
population = [{k: random.choice(v) if isinstance(v, list) else random.uniform(*v) for k,v in OPTIMIZATION_SPACE.items()} for _ in range(40)] |
|
|
best_overall_individual = None |
|
|
best_overall_fitness = -float('inf') |
|
|
for _ in range(80): |
|
|
fitnesses = [calculate_project_kpis(**ind)['irr'] for ind in population] |
|
|
if max(fitnesses) > best_overall_fitness: |
|
|
best_overall_fitness = max(fitnesses) |
|
|
best_overall_individual = population[np.argmax(fitnesses)] |
|
|
selected = [max(random.sample(list(zip(population, fitnesses)), 5), key=lambda i: i[1])[0] for _ in range(50)] |
|
|
next_gen = [] |
|
|
for i in range(0, 50, 2): |
|
|
p1, p2 = selected[i], selected[i+1] |
|
|
c1, c2 = copy.deepcopy(p1), copy.deepcopy(p2) |
|
|
if random.random() < 0.9: c1['technology'], c2['technology'] = p2['technology'], c1['technology'] |
|
|
if random.random() < 0.2: c1['export_market_mix'] = random.uniform(*OPTIMIZATION_SPACE['export_market_mix']) |
|
|
next_gen.extend([c1, c2]) |
|
|
population = next_gen |
|
|
kpis = calculate_project_kpis(**best_overall_individual) |
|
|
return {"Method": "Genetic Algorithm", **kpis, "Params": best_overall_individual} |
|
|
|
|
|
def run_bayesian_optimization(): |
|
|
print("Running Bayesian Optimization...") |
|
|
skopt_space = [ |
|
|
Real(OPTIMIZATION_SPACE['capacity_kta'][0], OPTIMIZATION_SPACE['capacity_kta'][1], name='capacity_kta'), |
|
|
Categorical(OPTIMIZATION_SPACE['technology'], name='technology'), |
|
|
Categorical(OPTIMIZATION_SPACE['sourcing_strategy'], name='sourcing_strategy'), |
|
|
Real(OPTIMIZATION_SPACE['export_market_mix'][0], OPTIMIZATION_SPACE['export_market_mix'][1], name='export_market_mix'), |
|
|
Categorical(OPTIMIZATION_SPACE['sell_byproducts'], name='sell_byproducts') |
|
|
] |
|
|
@use_named_args(skopt_space) |
|
|
def objective(**params): |
|
|
return -calculate_project_kpis(**params)['irr'] |
|
|
res = gp_minimize(objective, skopt_space, n_calls=90, random_state=42, n_initial_points=20) |
|
|
best_params = {space.name: val for space, val in zip(skopt_space, res.x)} |
|
|
kpis = calculate_project_kpis(**best_params) |
|
|
return {"Method": "Bayesian Opt", **kpis, "Params": best_params} |
|
|
|
|
|
def run_optuna_direct(): |
|
|
print("Running Optuna (TPE)...") |
|
|
def objective(trial): |
|
|
params = { |
|
|
"capacity_kta": trial.suggest_float("capacity_kta", *OPTIMIZATION_SPACE['capacity_kta']), |
|
|
"technology": trial.suggest_categorical("technology", OPTIMIZATION_SPACE['technology']), |
|
|
"sourcing_strategy": trial.suggest_categorical("sourcing_strategy", OPTIMIZATION_SPACE['sourcing_strategy']), |
|
|
"export_market_mix": trial.suggest_float("export_market_mix", *OPTIMIZATION_SPACE['export_market_mix']), |
|
|
"sell_byproducts": trial.suggest_categorical("sell_byproducts", OPTIMIZATION_SPACE['sell_byproducts']) |
|
|
} |
|
|
return calculate_project_kpis(**params)['irr'] |
|
|
study = optuna.create_study(direction="maximize") |
|
|
study.optimize(objective, n_trials=150, n_jobs=-1) |
|
|
kpis = calculate_project_kpis(**study.best_params) |
|
|
return {"Method": "Optuna (TPE - Direct)", **kpis, "Params": study.best_params} |
|
|
|
|
|
def display_and_save_results(df_results): |
|
|
print("\n--- Final Results and Comparison ---") |
|
|
df_display = pd.DataFrame() |
|
|
df_display['Method'] = df_results['Method'] |
|
|
df_display['Optimal IRR (%)'] = df_results['irr'].map('{:,.2f}%'.format) |
|
|
df_display['Annual Profit ($M)'] = (df_results['annual_profit'] / 1_000_000).map('{:,.1f}'.format) |
|
|
df_display['CAPEX ($M)'] = (df_results['total_capex'] / 1_000_000).map('{:,.1f}'.format) |
|
|
df_display['Payback (Yrs)'] = df_results['payback_period'].map('{:,.1f}'.format) |
|
|
param_df = pd.DataFrame(df_results['Params'].tolist()) |
|
|
param_df['capacity_kta'] = param_df['capacity_kta'].round(1) |
|
|
param_df['export_market_mix'] = (param_df['export_market_mix'] * 100).round(1).astype(str) + '%' |
|
|
df_display = pd.concat([df_display, param_df.rename(columns={ |
|
|
'capacity_kta': 'Capacity (KTA)', 'technology': 'Technology', 'sourcing_strategy': 'Sourcing', |
|
|
'export_market_mix': 'Export Mix', 'sell_byproducts': 'Sell Byproducts' |
|
|
})], axis=1) |
|
|
|
|
|
print("\n✅ **Final Comparison of Optimal Scenarios (Sorted by Best IRR)**") |
|
|
print("="*120) |
|
|
print(df_display.to_string(index=False)) |
|
|
print("="*120) |
|
|
df_display.to_csv("results.csv", index=False, encoding='utf-8-sig') |
|
|
|
|
|
def create_kpi_comparison_dashboard(df_results): |
|
|
print("\n--- Generating KPI Comparison Dashboard ---") |
|
|
df_plot = df_results.sort_values(by='irr', ascending=False) |
|
|
df_plot['annual_profit_M'] = df_plot['annual_profit'] / 1_000_000 |
|
|
df_plot['total_capex_M'] = df_plot['total_capex'] / 1_000_000 |
|
|
|
|
|
fig, axes = plt.subplots(2, 2, figsize=(20, 14)) |
|
|
fig.suptitle('Dashboard: Comprehensive Comparison of Optimization Methods', fontsize=22, weight='bold') |
|
|
palettes = ['viridis', 'plasma', 'magma', 'cividis'] |
|
|
metrics = [ |
|
|
('irr', 'Optimal IRR (%)', axes[0, 0]), |
|
|
('annual_profit_M', 'Annual Profit ($M)', axes[0, 1]), |
|
|
('total_capex_M', 'Total CAPEX ($M)', axes[1, 0]), |
|
|
('payback_period', 'Payback Period (Years)', axes[1, 1]) |
|
|
] |
|
|
|
|
|
for i, (metric, title, ax) in enumerate(metrics): |
|
|
sns.barplot(x=metric, y='Method', data=df_plot, ax=ax, palette=palettes[i]) |
|
|
ax.set_title(title, fontsize=16, weight='bold') |
|
|
ax.set_xlabel('') |
|
|
ax.set_ylabel('') |
|
|
for p in ax.patches: |
|
|
width = p.get_width() |
|
|
ax.text(width * 1.01, p.get_y() + p.get_height() / 2, |
|
|
f'{width:,.2f}', |
|
|
va='center', fontsize=11) |
|
|
|
|
|
plt.tight_layout(rect=[0, 0.03, 1, 0.95]) |
|
|
plt.savefig("static/images/kpi_dashboard.png", bbox_inches='tight') |
|
|
print("\n✅ A KPI dashboard graph has been saved as 'kpi_dashboard.png'") |
|
|
|
|
|
def run_sensitivity_analysis(best_params, base_irr): |
|
|
global PRODUCT_PRICES_USD_PER_TON, STRATEGY_DATA |
|
|
print("\n--- Sensitivity Analysis on Best Scenario ---") |
|
|
print(f"Analyzing sensitivity around the base IRR of {base_irr:,.2f}%") |
|
|
|
|
|
sensitivity_vars = { |
|
|
'Byproduct Caustic Soda Price': ('price', 'byproduct_caustic_soda'), |
|
|
'Sourcing Cost Integrated': ('strategy', 'Integrated_Production', 'sourcing_cost_per_ton_pvc'), |
|
|
'Domestic PVC S70 Price': ('price', 'pvc_s70_domestic') |
|
|
} |
|
|
variations = [-0.20, -0.10, 0.10, 0.20] |
|
|
results = [] |
|
|
original_prices = copy.deepcopy(PRODUCT_PRICES_USD_PER_TON) |
|
|
original_strategies = copy.deepcopy(STRATEGY_DATA) |
|
|
|
|
|
for key, path in sensitivity_vars.items(): |
|
|
for var in variations: |
|
|
PRODUCT_PRICES_USD_PER_TON = copy.deepcopy(original_prices) |
|
|
STRATEGY_DATA = copy.deepcopy(original_strategies) |
|
|
if path[0] == 'price': |
|
|
base_value = PRODUCT_PRICES_USD_PER_TON[path[1]] |
|
|
PRODUCT_PRICES_USD_PER_TON[path[1]] = base_value * (1 + var) |
|
|
else: |
|
|
base_value = STRATEGY_DATA[path[1]][path[2]] |
|
|
STRATEGY_DATA[path[1]][path[2]] = base_value * (1 + var) |
|
|
|
|
|
kpis = calculate_project_kpis(**best_params) |
|
|
results.append({ |
|
|
'Variable': key, 'Change': f'{var:+.0%}', |
|
|
'New IRR (%)': kpis['irr'], 'IRR Delta (%)': kpis['irr'] - base_irr |
|
|
}) |
|
|
|
|
|
PRODUCT_PRICES_USD_PER_TON = original_prices |
|
|
STRATEGY_DATA = original_strategies |
|
|
|
|
|
df_sensitivity = pd.DataFrame(results) |
|
|
print("\n✅ **Sensitivity Analysis Results**") |
|
|
print("="*60) |
|
|
print(df_sensitivity.to_string(index=False)) |
|
|
print("="*60) |
|
|
|
|
|
tornado_data = [] |
|
|
for var_name in df_sensitivity['Variable'].unique(): |
|
|
subset = df_sensitivity[df_sensitivity['Variable'] == var_name] |
|
|
min_delta = subset['IRR Delta (%)'].min() |
|
|
max_delta = subset['IRR Delta (%)'].max() |
|
|
tornado_data.append({ |
|
|
'Variable': var_name, |
|
|
'Min_Delta': min_delta, |
|
|
'Max_Delta': max_delta, |
|
|
'Range': max_delta - min_delta |
|
|
}) |
|
|
df_tornado = pd.DataFrame(tornado_data).sort_values('Range', ascending=True) |
|
|
|
|
|
fig, ax = plt.subplots(figsize=(12, 8)) |
|
|
y = np.arange(len(df_tornado)) |
|
|
ax.barh(y, df_tornado['Max_Delta'], color='mediumseagreen', label='Positive Impact') |
|
|
ax.barh(y, df_tornado['Min_Delta'], color='lightcoral', label='Negative Impact') |
|
|
ax.set_yticks(y) |
|
|
ax.set_yticklabels(df_tornado['Variable'], fontsize=12) |
|
|
ax.axvline(0, color='black', linewidth=0.8, linestyle='--') |
|
|
ax.set_title('Tornado Chart: IRR Sensitivity to Key Variables', fontsize=18, pad=20, weight='bold') |
|
|
ax.set_xlabel(f'Change in IRR (%) from Base IRR ({base_irr:.2f}%)', fontsize=14) |
|
|
ax.set_ylabel('Variable', fontsize=14) |
|
|
ax.legend() |
|
|
ax.grid(axis='x', linestyle='--', alpha=0.7) |
|
|
for i, (p, n) in enumerate(zip(df_tornado['Max_Delta'], df_tornado['Min_Delta'])): |
|
|
ax.text(p, i, f' +{p:.2f}%', va='center', ha='left', color='darkgreen') |
|
|
ax.text(n, i, f' {n:.2f}%', va='center', ha='right', color='darkred') |
|
|
plt.tight_layout() |
|
|
plt.savefig("static/images/sensitivity_analysis_tornado.png", bbox_inches='tight') |
|
|
print("\n✅ A sensitivity analysis Tornado chart has been saved as 'sensitivity_analysis_tornado.png'") |
|
|
|
|
|
if __name__ == "__main__": |
|
|
optimization_results = run_optimizations_without_ml() |
|
|
df_results = pd.DataFrame(optimization_results).sort_values(by="irr", ascending=False).reset_index(drop=True) |
|
|
df_results.round(2) |
|
|
display_and_save_results(df_results) |
|
|
create_kpi_comparison_dashboard(df_results) |
|
|
|
|
|
if not df_results.empty: |
|
|
best_scenario = df_results.iloc[0] |
|
|
run_sensitivity_analysis(best_scenario['Params'], best_scenario['irr']) |
