# Oil refinery production optimization (+PowerBI)¶

Description: In this document, we present an enhanced approach to oil refining optimization by integrating Power BI for improved decision-making and data visualization. For a full description of the model, you can read more about it here.

This notebook showcases how Power BI, when combined with mathematical optimization, delivers powerful, actionable insights. By leveraging the dynamic reporting features of Power BI and the optimization capabilities of AMPL, users can interact with data in real-time and make well-informed decisions more effectively.

Tags: Oil production, Production optimization, Profitability, Refinery, mip, highs, powerbi, industry, scheduling, data-science, data-analysis, decision-making

Notebook author: Mikhail Riabtsev <mail@solverytic.com>

## 1. Introduction¶

In this document, we present an enhanced approach to oil refining optimization by integrating Power BI for improved decision-making and data visualization. For a full description of the model, you can read more about it here.

This notebook showcases how Power BI, when combined with mathematical optimization, delivers powerful, actionable insights. By leveraging the dynamic reporting features of Power BI and the optimization capabilities of AMPL, users can interact with data in real-time and make well-informed decisions more effectively. This innovative method offers several key benefits:

**Interactive Dashboards:**Power BI’s user-friendly interface allows easy exploration of scenarios and results, without requiring advanced technical knowledge.**Real-time Data Integration**: Dynamic data inputs ensure optimization is always based on the most current information, keeping decisions timely and relevant.**Scalability**: Power BI’s ability to handle large datasets makes it ideal for implementing optimization models on a larger scale.**User-Friendly Reporting**: Clear and intuitive visual representations of optimization results enhance understanding, facilitating smoother communication across teams. Together, AMPL and Power BI create a streamlined process, bridging the gap between model development and real-world decision-making.

## 1. Download Necessary Extensions and Libraries¶

Ensure you have a registered Power BI Service account and uploaded Oil_refining.pbix file to your Power BI Service.

```
# Install dependencies
%pip install -q amplpy powerbiclient
from powerbiclient import Report, models
from io import StringIO
from ipywidgets import interact
import requests
import pandas as pd
import numpy as np
#import matplotlib.pyplot as plt
```

```
# Google Colab & Kaggle integration
from amplpy import AMPL, ampl_notebook
ampl = ampl_notebook(
modules=["highs"], # modules to install
license_uuid="default", # license to use
) # instantiate AMPL object and register magics
```

## 3. Embed the Power BI Report with Data¶

Integrate your Power BI report to Jupyter Notebook.

```
# Get the group and report ID from PowerBI service
group_id="c498f9cc-fb27-4ba4-96b6-d4c1813624c2" # YOU HAVE TO PUT HERE YOUR POWER BI GROUP ID OR WORKSPACE ID
report_id="7500c9bc-ac5e-4e91-a8ff-bba318f3d6a9" # YOU HAVE TO PUT HERE YOUR POWER BI REPORT ID
# Creating an instance of the Report class by passing group_id, report_id, and interactive_auth
report = Report(group_id=group_id, report_id=report_id, auth = interactive_auth)
#report.set_size(800, 1240) # This line is commented out but would set the size of the report's output.
report # Show the power BI report with data.
```

## 4. AMPL Model Formulation¶

Define and set up the optimization model.

```
%%ampl_eval
reset;
### SETS
set CRUD ; # Types of crude oil (Crude1, Crude2)
set DIST ; # Distillation products
set REF ; # Reforming products
set CRACK ; # Cracking products
set LUBR ; # Lubricating products
set PROD ; # Final products
set STAT ; # Delivery stations
set D_MODE ; # Operating modes of the equipment of the Distillation installation
set R_MODE ; # Operating modes of the equipment of the Reforming installation
set CR_MODE ; # Operating modes of the equipment of the Cracking installation
set L_MODE ; # Operating modes of the equipment of the Lubrication process
set POLLUT ; # Types of pollutants
set DISTILLATION within {D_MODE, CRUD, DIST} ; # Pairs of the Distillation process
set REFORMING within {R_MODE, DIST, REF} ; # Pairs of the Reforming process
set DIST_R := setof{(m,d,r) in REFORMING}(d) ; # List of included components
set CRACKING within {CR_MODE, DIST, CRACK} ; # Pairs of the Cracking process
set DIST_CR := setof{(m,d,cr) in CRACKING}(d) ; # List of included components
set LUBRICATING within {L_MODE, DIST, LUBR} ; # Pairs of the Lubricating process
set DIST_L := setof{(m,d,l) in LUBRICATING}(d) ; # List of included components
set BLENDING within # Pairs of Intermediate and final products involved in blending
{DIST union REF union CRACK union LUBR, PROD} ;
set INTERMED:= setof{(i,j) in BLENDING}i ; # Set of Intermediate products before blending
### PARAMETERS
param nPeriod >= 0 ; # Number of weeks in the planning period
param nPeriodByYear >= 0 ; # Number of nPeriods in the Year
## Crude oil
param crude_Min{CRUD} = 0 ; # Minimum supply limits for each crude oil type
param crude_Max_capacity{c in CRUD} >= crude_Min[c] ; # Maximum supply limits for each crude oil type
param crude_Cost{CRUD, 1..nPeriod} >= 0 ; # Cost of crude oil per week
## Distillation
param distill_Yield{DISTILLATION} >= 0 ; # Yield of products
param distill_Pollute{D_MODE, CRUD, POLLUT} ; # Pollutant emissions
param distill_Cost{D_MODE, CRUD} >= 0 ; # Cost of process
param distill_Waste_Cost{D_MODE, CRUD} >= 0 ; # Residue disposal cost
param distill_Equipment_Setup_Period{D_MODE} >= 0 ; # Equipment setup period
param distill_Equipment_Setup_Cost{D_MODE} >= 0 ; # Equipment setup cost
param distill_Max_capacity{D_MODE} >= 0 ; # Maximum capacity
## Reforming
param reform_Yield{REFORMING} >= 0 ; # Yield of products
param reform_Pollute{R_MODE, DIST_R, POLLUT} ; # Pollutant emissions
param reform_Cost{R_MODE, DIST_R} >= 0 ; # Cost of process
param reform_Waste_Cost{R_MODE, DIST_R} >= 0 ; # Residue disposal cost
param reform_Max_capacity{R_MODE} >= 0 ; # Maximum capacity
param reform_Equipment_Setup_Period{R_MODE} >= 0 ; # Equipment setup period
param reform_Equipment_Setup_Cost{R_MODE} >= 0 ;# Equipment setup cost
## Cracking
param crack_Yield{CRACKING} >= 0 ; # Yield of products
param crack_Pollute{CR_MODE, DIST_CR, POLLUT} ; # Pollutant emissions
param crack_Cost{CR_MODE, DIST_CR} >= 0 ; # Cost of process
param crack_Waste_Cost{CR_MODE, DIST_CR} >= 0 ; # Residue disposal cost
param crack_Max_capacity{CR_MODE} >= 0 ; # Maximum capacity
param crack_Equipment_Setup_Period{CR_MODE} >= 0 ; # Equipment setup period
param crack_Equipment_Setup_Cost{CR_MODE} >=0 ; # Equipment setup cost
## Lubrication
param lube_Yield{LUBRICATING} >= 0 ; # Yield of products
param lube_Pollute{L_MODE, DIST_L, POLLUT} ; # Pollutant emissions
param lube_Waste_Cost{L_MODE, DIST_L} >= 0 ; # Residue disposal cost
param lube_Cost{L_MODE, DIST_L} >= 0 ; # Cost of process
param lube_Max_capacity{L_MODE} >= 0 ; # Maximum production of lube oil
param lube_Equipment_Setup_Period{L_MODE} >= 0 ;# Equipment setup period
param lube_Equipment_Setup_Cost{L_MODE} >= 0 ; # Equipment setup cost
param lube_limit_Min = 0 ; # Minimum production of lube oil
param lube_limit_Max >= lube_limit_Min ; # Maximum production of lube oil
## Intermediate components
param Intermed_Octane{INTERMED} >= 0 ; # Octane number
param Intermed_VaporPressure{INTERMED} >= 0 ; # Vapor pressure
## Blending
param blending_Cost{PROD} >= 0 ; # Cost of blending
## Products
param prod_Octane_Min{PROD} >= 0 ; # Minimum octane number for final products
param prod_VaporPressure_Max{PROD} >= 0 ; # Maximum vapor pressure for final products
param prod_Premium_Regular_Gas_Min >= 0 ; # Minimum production of premium gas
param prod_FuelOil_Ratio{INTERMED} >= 0 ; # Ratios for fuel oil production components
## Storage
param storage_Capacity{PROD} >=0 ; # Storage capacity for each product
param storage_Cost{PROD} >= 0 ; # Storage cost per product
param storage_Waste{PROD} >= 0 ; # Waste during storage
## Product delivery
param delivery_Cost{PROD, STAT} >= 0 ; # Delivery cost per product to each station
## Plant
param plant_Shutdown_Period >= 0 ; # Equipment setup period
param plant_Shutdown_Cost >= 0 ; # Equipment setup cost
param plant_Const_Cost >= 0 ; # Plant fixed costs
## Market
param seasonal_Base_Demand{PROD, 1..nPeriod} >= 0 ; # Base demand for products per week.
param seasonal_Base_Price {PROD, 1..nPeriod} >= 0 ; # Base price for products per week.
# Price elasticity
param nStep integer > 0 ; # Number of steps for price elasticity
param price_nStep_Value{1..nStep+1} ; # Step values for price elasticity
param demand_nStep_Value{1..nStep+1} ; # Step values for price elasticity
## Finance
param discount_Rate >= 0 ; # Discount rate for future cash flows
param initial_Cash >= 0 ; # Initial cash available
## Loans
set LOANS; # Set of loan periods
set LOAN_param; # Parameters of loans (term, interest, Max_Money)
param loan{LOANS, LOAN_param} >= 0 ; # Conditions for obtaining credit
### VARIABLES
## Plant working
var Plant_Working{t in 1..nPeriod} binary; # 1 if the plant is running. 0 if the plant is shutdown
## Distillation
var Crude_Supply{D_MODE, CRUD, 1..nPeriod} >= 0 ; # Amount of crude supplied
var Distill_X{D_MODE, 1..nPeriod} binary ; # Additional binary variable for selecting the operating mode of the equipment
## Reforming
var Distill_to_Reforming{R_MODE, DIST_R, 1..nPeriod} >= 0 ; # Quantity of distillation products used for Reforming
var Reform_X{R_MODE, 1..nPeriod} binary ; # Additional binary variable for selecting the operating mode of the equipment
## Cracking
var Distill_to_Cracking{CR_MODE, DIST_CR, 1..nPeriod} >= 0 ;# Quantity of distillation products used for Cracking
var Cracking_X{CR_MODE, 1..nPeriod} binary ; # Additional binary variable for selecting the operating mode of the equipment
## Lubricating
var Distill_to_Lubricating{L_MODE, DIST_L, 1..nPeriod} >= 0 ;# Quantity of distillation products used for Lubricating
var Lubricating_X{L_MODE, 1..nPeriod} binary ; # Additional binary variable for selecting the operating mode of the equipment
## Blending:
var Blending{BLENDING, 1..nPeriod} >= 0 ; # Amount of ingredients mixed to obtain final products
## Demand
var Demand{PROD, STAT, 1..nPeriod, 1..nStep} >= 0 ; # Demand for each pr oduct at each station over time
var X{PROD, STAT, 1..nPeriod, 1..nStep} binary ; # Additional binary variable for demand steps (1 if for product p in period t the price is selected at step nStep, or 0 otherwise)
## Storage
var Storage_Fraction{p in PROD, t in 1..nPeriod} = # Amount of each product in Storage each period
sum{tt in 1..t} (sum{(i,p) in BLENDING} Blending[i,p,tt]
- sum{s in STAT, n in 1..nStep} Demand[p,s,tt,n]) * (1-storage_Waste[p]/100) ;
## Loan
# Amount of loan taken every period
var Loan_In{l in LOANS, 1..nPeriod-1} >= 0 , <= loan[l,'Max_Money'] ;
# Amount of loan taken every period
var Loan_Out{l in LOANS, t in 2..nPeriod} = Loan_In[l,t-1] * (1+(loan[l, 'interest']/nPeriodByYear)/100);
## Pollutant emissions
var Waste_Pollutant{p in POLLUT, t in 1..nPeriod} =
# Pollution from waste disposal from the Distillation process
sum{m in D_MODE, c in CRUD} Crude_Supply[m,c,t] * (1 - sum{d in DIST}distill_Yield[m,c,d]) * distill_Pollute[m,c,p]
# Pollution from waste disposal from the Reforming process
+ sum{m in R_MODE, c in DIST_R} Distill_to_Reforming[m,c,t] * (1 - sum{d in REF} reform_Yield[m,c,d]) * reform_Pollute[m,c,p]
# Pollution from waste disposal from the Cracking process
+ sum{m in CR_MODE, c in DIST_CR} Distill_to_Cracking[m,c,t] * (1 - sum{d in CRACK} crack_Yield[m,c,d]) * crack_Pollute[m,c,p]
# Pollution from waste disposal from the Lubricating process
+ sum{m in L_MODE, c in DIST_L} Distill_to_Lubricating[m,c,t] * (1 - sum{d in LUBR} lube_Yield[m,c,d]) * lube_Pollute[m,c,p];
## Cash flow with incomes and costs
var CashFlow{t in 1..nPeriod} =
# sales income
sum{p in PROD, s in STAT, n in 1..nStep} Demand[p,s,t,n] * (seasonal_Base_Price[p,t] * (1 + price_nStep_Value[n]/100) - delivery_Cost[p,s])
# minus the cost of purchasing crude oil + the costs of the Distillation process + costs for disposal of waste from the Distillation process
- sum{m in D_MODE, c in CRUD} Crude_Supply[m,c,t] * (crude_Cost[c,t] + distill_Cost[m,c] + (1 - sum{d in DIST}distill_Yield[m,c,d]) * distill_Waste_Cost[m,c])
# minus the costs of the Reforming process + costs for disposal of waste from the Reforming process
- sum{m in R_MODE, c in DIST_R} Distill_to_Reforming[m,c,t] * (reform_Cost[m,c] + (1 - sum{d in REF} reform_Yield[m,c,d]) * reform_Waste_Cost[m,c])
# minus the costs of the Cracking process + costs for disposal of waste from the Cracking process
- sum{m in CR_MODE, c in DIST_CR} Distill_to_Cracking[m,c,t] * (crack_Cost[m,c] + (1 - sum{d in CRACK} crack_Yield[m,c,d]) * crack_Waste_Cost[m,c])
# minus the costs of the Lubricating process + costs for disposal of waste from the Lubricating process
- sum{m in L_MODE, c in DIST_L} Distill_to_Lubricating[m,c,t] * (lube_Cost[m,c] + (1 - sum{d in LUBR} lube_Yield[m,c,d]) * lube_Waste_Cost[m,c])
# minus the costs of reconfiguring equipment (if available)
- (if t > 1 then /*Nonliner piece*/
# Distillation equipment
sum{m in D_MODE} (if Distill_X[m,t] - Distill_X[m,t-1] > 0 then distill_Equipment_Setup_Cost[m] else 0)
# Reforming equipment
+ sum{m in R_MODE} (if Reform_X[m,t] - Reform_X[m,t-1] > 0 then reform_Equipment_Setup_Cost[m] else 0)
# Cracking equipment
+ sum{m in CR_MODE} (if Cracking_X[m,t] - Cracking_X[m,t-1] > 0 then crack_Equipment_Setup_Cost[m] else 0)
# Lubricating equipment
+ sum{m in L_MODE}(if Lubricating_X[m,t] - Lubricating_X[m,t-1] > 0 then lube_Equipment_Setup_Cost[m] else 0)
else 0)
# minus the cost of the Blending
- sum{(i,p) in BLENDING} Blending[i,p,t] * blending_Cost[p]
# minus cost of shutdown
- (1 - Plant_Working[t]) * plant_Shutdown_Cost
# minus the storage cost
- sum{p in PROD} Storage_Fraction[p,t] * storage_Cost[p]
# minus fixed plant costs
- plant_Const_Cost
# plus the amount of the loan received
+ (if t = nPeriod then 0 else sum{l in LOANS} Loan_In[l,t])
# minus the amount of repaid loans with accrued interest
- (if t = 1 then 0 else sum{l in LOANS} Loan_Out[l,t]);
### OBJECTIVE FUNCTION
# Maximize the total profit considering all incomes and costs, discounted by the discount rate
maximize Total_Profit:
sum{t in 1..nPeriod} CashFlow[t] / (1 + (discount_Rate/nPeriodByYear)/100)^t ;
### CONSTRAINTS
subject to
## Distillation
# Ensure that total supply of crude oil does not exceed the maximum capacity.
Crude_Supply_Min_Max{c in CRUD, t in 1..nPeriod}:crude_Min[c] <= sum{m in D_MODE} Crude_Supply[m,c,t] <= crude_Max_capacity[c];
# Ensure distillation capacity does not exceed the maximum. Сapacity is reduced during downtime caused by equipment reconfiguring.
DistillCapacity_Max {m in D_MODE, t in 1..nPeriod}: sum{c in CRUD} Crude_Supply[m,c,t] <= distill_Max_capacity [m]
/* Nonliner piece*/ * (if t > 1 then (if Distill_X[m,t] > Distill_X[m,t-1] then 1 - distill_Equipment_Setup_Period[m] else 1) else 1) ;
# Ensure only one mode per period
ForEachPeriodOnlyOne_Distill_X { t in 1..nPeriod}: sum{m in D_MODE} Distill_X[m,t] <= 1 ;
# Ensure that Crude_Supply > 0 only when Distill_X > 0
Distill_{m in D_MODE, t in 1..nPeriod}: sum{c in CRUD} Crude_Supply[m,c,t] <= Distill_X[m,t] * 10e5;
## Reforming
# Ensure reforming capacity does not exceed the maximum. Сapacity is reduced during downtime caused by equipment reconfiguring.
Reforming_Capacity_Max {m in R_MODE, t in 1..nPeriod}: sum{c in DIST_R} Distill_to_Reforming[m,c,t] <= reform_Max_capacity[m]
/* Nonliner piece*/ * (if t > 1 then (if Reform_X[m,t] > Reform_X[m,t-1] then 1 - reform_Equipment_Setup_Period[m] else 1) else 1);
# Ensure only one mode per period
ForEachPeriodOnlyOne_Reform_X { t in 1..nPeriod}: sum{m in R_MODE} Reform_X[m,t] <= 1 ;
# Ensure that Distill_to_Reforming > 0 only when Reform_X > 0
Reform_{m in R_MODE, t in 1..nPeriod}:sum{c in DIST_R}Distill_to_Reforming[m,c,t] <= Reform_X[m,t] * 10e5;
## Cracking
# Ensure cracking capacity does not exceed the maximum. Сapacity is reduced during downtime caused by equipment reconfiguring.
CrackingCapacity_Max {m in CR_MODE, t in 1..nPeriod}: sum{c in DIST_CR} Distill_to_Cracking[m,c,t] <= crack_Max_capacity[m]
/* Nonliner piece*/ * (if t > 1 then (if Cracking_X[m,t] > Cracking_X[m,t-1] then 1 - crack_Equipment_Setup_Period[m] else 1) else 1);
# Ensure only one mode per period
ForEachPeriodOnlyOne_Cracking_X { t in 1..nPeriod}: sum{m in CR_MODE} Cracking_X[m,t] <= 1 ;
# Ensure that Distill_to_Cracking > 0 only when Cracking_X > 0
Cracking_{m in CR_MODE, t in 1..nPeriod}: sum{c in DIST_CR} Distill_to_Cracking[m,c,t] <= Cracking_X[m,t] * 10e5;
## Lubricating
# Ensure cracking does not exceed the minimum & maximum volume
Lube_Oil_Min_Max {t in 1..nPeriod}:
lube_limit_Min <= sum{(m,d,l)in LUBRICATING} Distill_to_Lubricating[m,d,t] * lube_Yield[m,d,l] <= lube_limit_Max ;
# Сapacity is reduced during downtime caused by equipment reconfiguring.
Lube_Oil_Capacity_Max {m in L_MODE, t in 1..nPeriod}: sum{c in DIST_L} Distill_to_Lubricating[m,c,t] <= lube_Max_capacity[m]
/* Nonliner piece*/ * (if t > 1 then (if Lubricating_X[m,t] > Lubricating_X[m,t-1] then 1 - lube_Equipment_Setup_Period[m] else 1) else 1); # Nonliner piece
# Ensure only one mode per period
ForEachPeriodOnlyOne_Lubricating_X { t in 1..nPeriod}: sum{m in L_MODE} Lubricating_X[m,t] <= 1 ;
# Ensure that Distill_to_Lubricating> 0 only when Lubricating_X > 0
Lubricating_{m in L_MODE, t in 1..nPeriod}: sum{c in DIST_L}Distill_to_Lubricating[m,c,t] <= Lubricating_X[m,t] * 10e5;
## Blending
PremiumRegularGasRatio {t in 1..nPeriod}: # Premium gasoline production: at least 40% of regular gasoline production
sum{(i,p) in BLENDING: p='Premium Gasoline'} Blending[i,p,t] >=
sum{(i,p) in BLENDING: p='Regular Gasoline'} Blending[i,p,t] * prod_Premium_Regular_Gas_Min/100 ;
OctaneNumberMin {p in PROD, t in 1..nPeriod: p = 'Premium Gasoline'}:# Ensure octane number requirements for final products
sum{(i,p) in BLENDING} Blending[i,p,t] * Intermed_Octane[i] >=
sum{(ii,p) in BLENDING} Blending[ii,p,t] * prod_Octane_Min[p] ;
OctaneNumberMin_2 {p in PROD, t in 1..nPeriod: p = 'Regular Gasoline'}:# Ensure octane number requirements for final products
sum{(i,p) in BLENDING} Blending[i,p,t] * Intermed_Octane[i] >=
sum{(ii,p) in BLENDING} Blending[ii,p,t] * prod_Octane_Min[p] ;
VaporPressure_Max {p in PROD, t in 1..nPeriod: p='Jet Fuel'}:# Ensure vapor pressure limits for products
sum{(i,p) in BLENDING} Blending[i,p,t] * Intermed_VaporPressure[i] <=
sum{(i,p) in BLENDING} Blending[i,p,t] * prod_VaporPressure_Max[p] ;
FuelOilRatio {(i,p) in BLENDING, t in 1..nPeriod: p = 'Fuel Oil'}: # Maintain the correct ratio for fuel oil production
Blending['Residuum',p, t] = Blending[i,p,t] / prod_FuelOil_Ratio[i] ;
## Storage
# Ensure storage fractions are non-negative and within capacity
Storage_non_negative {p in PROD, t in 1..nPeriod}:
0 <= Storage_Fraction[p,t] <= storage_Capacity[p];
## Financial calculations
CashFlow_Balance{t in 1..nPeriod}: # Cash flow balance
initial_Cash
+ sum{tt in 1..t}CashFlow[tt] >= 0;
## Ensure sufficient distillation output
Distillation_Out {d in DIST, t in 1..nPeriod}:
sum{m in D_MODE, c in CRUD} Crude_Supply[m,c,t] * distill_Yield[m,c,d]
>=
sum{m in R_MODE, dd in DIST_R: dd=d} Distill_to_Reforming[m,dd,t]
+ sum{m in CR_MODE, dd in DIST_CR: dd=d} Distill_to_Cracking[m,dd,t]
+ sum{m in L_MODE, dd in DIST_L: dd=d} Distill_to_Lubricating[m,dd,t] ;
## Balance INTERMED products before blending
INTERMED_Balance{i in INTERMED, t in 1..nPeriod}:
sum{(i,p) in BLENDING} Blending[i,p,t] <=
sum{(m,c,d) in DISTILLATION:d=i} Crude_Supply[m,c,t] * distill_Yield[m,c,d]
- sum{m in R_MODE, d in DIST_R: d=i} Distill_to_Reforming[m,d,t]
- sum{m in CR_MODE, d in DIST_CR: d=i} Distill_to_Cracking[m,d,t]
- sum{m in L_MODE, d in DIST_L: d=i} Distill_to_Lubricating[m,d,t]
+ sum{(m,d,r) in REFORMING: r=i} Distill_to_Reforming[m,d,t] * reform_Yield[m,d,r]
+ sum{(m,d,cr) in CRACKING: cr=i} Distill_to_Cracking[m,d,t] * crack_Yield[m,d,cr]
+ sum{(m,d,l) in LUBRICATING: l=i} Distill_to_Lubricating[m,d,t] * lube_Yield[m,d,l] ;
## Prices elastisity
ForEachPeriodOnlyOne_X {p in PROD, s in STAT, t in 1..nPeriod}:
sum{i in 1..nStep} X [p,s,t,i] = 1 ; # Ensure only one price per product per period
# Double restriction. Ensure demand is within the range dictated by price steps.
# The constraint of type min must be set after solving a model without this constraint. #Use Min constraint after ensuring that it does not conflict with other model constraints.
Demand_Min {p in PROD, s in STAT, t in 1..nPeriod, n in 1..nStep: p!='Lube Oil'}:
Demand[p,s,t,n] >= X[p,s,t,n] * seasonal_Base_Demand[p,t] * (1+demand_nStep_Value[n]/100) ;
Demand_Max {p in PROD, s in STAT, t in 1..nPeriod, n in 1..nStep}:
Demand[p,s,t,n] <= X[p,s,t,n] * seasonal_Base_Demand[p,t] * (1+demand_nStep_Value[n+1]/100) ;
## Shutdown contstraints
# The number of working periods is equal to the total number of periods under consideration minus the duration of the planned suspension
Plant_Working_nPeriods: sum{t in 1..nPeriod - plant_Shutdown_Period + 1} Plant_Working[t] = nPeriod - plant_Shutdown_Period;
# Additionally, there is a restriction in case of a longer (more than 1 period) plant shutdown.
Shutdown_Distill{t in 1..nPeriod - plant_Shutdown_Period + 1}:
sum{m in D_MODE, c in CRUD, tt in 0..plant_Shutdown_Period-1} Crude_Supply[m,c,t+tt] <= Plant_Working[t] * 10e5;
```

## 5. Data change¶

You can change model parameters directly in a Power BI report.

## 5. Download Data from Power BI Report¶

Retrieve the required data for analysis.

```
pages = report.get_pages() # Retrieve the list of pages from the report object, storing them in the 'pages' variable.
# The following lines are commented out but would iterate through the list of pages if uncommented.
#for page in range(len(pages)): # Loop over each page index from 0 to the total number of pages.
#print("displayName:", pages[page]['displayName']," name:", pages[page]['name']) # For each page, print its 'displayName' (the human-readable name) and 'name' (the internal identifier).
range_of_visual = { # List visuals for import data to AMPL parameters
'a6ee369779ca1dbbe625': ['CRUD','crude_Max_capacity', 'crude_Cost'],
'ReportSectione9e707e99075fef5f637': ['DIST', 'POLLUT', 'D_MODE', 'distill_Waste_Cost', 'distill_Max_capacity','distill_Equipment_Setup_Period','distill_Pollute','distill_Cost', 'distill_Yield', 'distill_Equipment_Setup_Cost'],
'1f30731c20476b22e3b8': ['REF', 'R_MODE', 'reform_Max_capacity', 'reform_Equipment_Setup_Period', 'reform_Equipment_Setup_Cost', 'reform_Yield', 'reform_Cost', 'reform_Pollute', 'reform_Waste_Cost'],
'baaf6d948d8009021852': ['CRACK', 'CR_MODE', 'crack_Yield', 'crack_Pollute', 'crack_Cost', 'crack_Waste_Cost', 'crack_Max_capacity', 'crack_Equipment_Setup_Period', 'crack_Equipment_Setup_Cost'],
'df4105df5589821b0042': ['LUBR', 'L_MODE', 'lube_Max_capacity', 'lube_Equipment_Setup_Period', 'lube_Equipment_Setup_Cost', 'lube_limit_Max', 'lube_Yield', 'lube_Pollute', 'lube_Cost', 'lube_Waste_Cost'],
'5f6474c4e0a5741fa523': ['PROD', 'prod_FuelOil_Ratio', 'prod_VaporPressure_Max', 'prod_Premium_Regular_Gas_Min', 'blending_recipe', 'prod_Octane_Min', 'Intermed_VaporPressure', 'Intermed_Octane', 'blending_Cost'],
'824aaabf4036f0ccfbb2': ['STAT', 'storage_Capacity', 'storage_Cost', 'storage_Waste', 'delivery_Cost'],
'a0c8f0fc396f1f12c937': ['plant_Const_Cost', 'plant_Shutdown_Period', 'plant_Shutdown_Cost', 'initial_Cash', 'discount_Rate'],
'015cea0c3cc0be908cc3': ['LOAN_param','LOANS','loan'],
'eac0f7120065aad83ca0': ['demand_nStep_Value', 'seasonal_Base_Price', 'seasonal_Base_Demand', 'price_nStep_Value'],
}
ampl_ = dict() # Initialize an empty dictionary to store data for each visual param element
for key_s in range_of_visual: # Loop through each page in the list of Power BI pages
report.set_active_page(key_s) # Set the current page as active in the Power BI report
visuals = report.visuals_on_page(key_s) # Get the list of visual elements on the active page
for i in range(len(range_of_visual[key_s])): # Loop through each visual element on the active page
#print(range_of_visual[key_s][i]) # Debugging statement to print the current visual element name
visual = next(filter(lambda visual: visual['title'] == range_of_visual[key_s][i], visuals)) # Find the visual element from the active page's visuals list that matches the name in range_of_visual_param
summarized_exported_data = report.export_visual_data(key_s, visual['name'], rows=1000, export_data_type=models.ExportDataType.SUMMARIZED.value) # Export summarized data from the visual element
data = StringIO(summarized_exported_data) # Convert the exported data to a string format for processing
ampl_[range_of_visual[key_s][i]] = pd.read_csv(data) # Read the data as a CSV and store it in the dictionary
# After looping through all pages and visuals, reset to the main page
report.set_active_page(pages[0]['name']) # Re-activate the main page
```

## 6. Load Data into AMPL¶

Seamlessly import your data into the AMPL system.

```
### SETS
ampl.set['DISTILLATION'] = ampl_['distill_Yield'].apply(lambda row: (row['D_MODE'], row['CRUD'], row['DIST']), axis=1)
ampl.set['REFORMING'] = ampl_['reform_Yield'].apply(lambda row: (row['R_MODE'], row['DIST_R'], row['REF']), axis=1)
ampl.set['CRACKING'] = ampl_['crack_Yield'].apply(lambda row: (row['CR_MODE'], row['DIST_CR'], row['CRACK']), axis=1)
ampl.set['LUBRICATING'] = ampl_['lube_Yield'].apply(lambda row: (row['L_MODE'], row['DIST_L'], row['LUBR']), axis=1)
ampl.set['BLENDING'] = ampl_['blending_recipe'].apply(lambda row: (row['INTERMED'], row['PROD']), axis=1)
### PARAMETERS
ampl.param['nPeriod'] = len(set(ampl_['crude_Cost']['nPeriod']))
ampl.param['nStep'] = len(set(ampl_['demand_nStep_Value']['nStep']))-1
ampl.param['nPeriodByYear'] = 12
sets_ = ['CRUD','DIST', 'POLLUT', 'D_MODE', 'REF', 'R_MODE', 'CRACK', 'CR_MODE', 'LUBR', 'L_MODE', 'PROD', 'LOAN_param','LOANS', 'STAT']
for key_dict in ampl_.keys(): # Loop through all keys in the 'ampl_' dictionary
if key_dict in sets_: # Check if the current key is one of the predefined sets
if len(ampl_[key_dict]) > 1: # If the DataFrame with more than one row
ampl.set[key_dict] = set(ampl_[key_dict].squeeze()) # Set the parameter in 'ampl'
else: # If only one row, directly extract the value and create a set
ampl.set[key_dict] = set(ampl_[key_dict][key_dict])
else:
df = ampl_[key_dict] # Get the dataframe corresponding to the current key from the 'ampl_param' dictionary
set_of_columns = df.loc[:, df.columns != key_dict].columns.tolist() # Get a list of indexing columns (excluding column 'key_var')
if set_of_columns: # If there are any remaining columns in 'r_df'
if key_dict not in {'blending_recipe'}: # Check if the current key is not 'blending_recipe'. we use 'blending_recipe' only for set
ampl.param[key_dict] = df.set_index(set_of_columns)[key_dict] # Set the parameter in 'ampl' with 'key_var' as the index based on the remaining columns
else:
ampl.param[key_dict] = ampl_[key_dict].squeeze() # Assign the squeezed version (i.e., reduced to lowest dimension) of 'ampl_param[key_var]' to 'ampl.param'
```

## 7. Solve the Problem¶

Run the optimization model to get your solution

```
%%ampl_eval
option solver highs; # Choosing a solver
# Defining Output Settings
option show_stats 1; # (1) Show statistical information about the size of the problem. Default 0 (statistics are not displayed)
option display_1col 0; # Data Display Settings
option omit_zero_rows 1; # Hide rows with 0 values. Default (0)
solve; # Solve the model
```

## 6. Display the solution¶

```
%%ampl_eval
# Display the values of the variables:
display Plant_Working, Crude_Supply, Distill_to_Reforming, Distill_to_Cracking, Distill_to_Lubricating, Blending, OctaneNumberMin, VaporPressure_Max, Storage_Fraction, Demand, Loan_In, Loan_Out, CashFlow, Waste_Pollutant;
# Calculate and print the octane number for each product over all periods.
printf {p in PROD, t in 1..nPeriod}: "Octane number for products: %u %s %6.2f\n", t, p, sum{(i,p) in BLENDING} Blending[i,p,t] * Intermed_Octane[i] / sum{(ii,p) in BLENDING}(if Blending[ii,p,t] = 0 then 1 else Blending[ii,p,t]);
# Calculate and print the vapor pressure limits for "Jet Fuel" over all periods.
printf {p in PROD, t in 1..nPeriod: p='Jet Fuel'}: "Vapor pressure limits of product Jet Fuel: %u %s %6.2f\n", t, p, sum{(i,p) in BLENDING} Blending[i,p,t] * Intermed_VaporPressure[i] / (if sum{(ii,p) in BLENDING} Blending[ii,p,t] = 0 then 1 else sum{(ii,p) in BLENDING}Blending[ii,p,t]) ;
# Calculate and print the percentage ratio of "Premium Gasoline" to "Regular Gasoline" over all periods.
printf {t in 1..nPeriod}: "Premium / Regular percentage: %u %6.2f\n", t,
sum{(i,p) in BLENDING: p ='Premium Gasoline'} Blending[i,p,t] / (if sum{(ii,pp) in BLENDING: pp ='Regular Gasoline'}Blending[ii,pp,t] = 0 then 1 else sum{(ii,pp) in BLENDING: pp ='Regular Gasoline'}Blending[ii,pp,t]) ;
# Calculate and print the production ratio for "Fuel Oil" for each blending component over all periods.
printf {(i,p) in BLENDING, t in 1..nPeriod: p='Fuel Oil'}: "Production Fuel Oil ratio: %u %s %6.2f\n", t, i,
Blending[i,p,t] / (if Blending['Residuum',p,t] = 0 then 1 else Blending['Residuum',p,t]);
```

## 7. Retrieve solution and export it to a *.xlsx file¶

Write the data to an *.xlsx file.

```
! pip install openpyxl
amplvar = dict()
# Create an ExcelWriter object for saving DataFrames to an Excel file at the specified path
with pd.ExcelWriter('oil_refyning.xlsx') as writer:
# Generate a list of all variable names from the AMPL model
list_of_ampl_variables = [item[0] for item in ampl.get_variables()]
# Iterate over each variable name in the list
for key_ampl in list_of_ampl_variables:
# Skip certain variables that are not to be processed (these variables won't be included in the output)
if key_ampl not in ['X', 'Cracking_X', 'Distill_X', 'Lubricating_X', 'Reform_X']:
# Convert the AMPL variable data to a pandas DataFrame
df = ampl.var[key_ampl].to_pandas()
# Filter the DataFrame to include only rows where the variable's value is greater than a small threshold (10^-5)
filtered_df = df[df[f"{key_ampl}.val"] > 10e-5]
# Save the filtered DataFrame to the corresponding sheet in the Excel file
if not filtered_df.empty: # Ensure that only non-empty DataFrames are saved
filtered_df.to_excel(writer, sheet_name=key_ampl, index=True)
# Save the filtered DataFrame to a dictionary
#amplvar[key_ampl] = filtered_df
```

## 8. Retrieve solution and export it to a Google Sheets document¶

Write the data to an Google Sheets document.

```
! pip install google-auth google-auth-oauthlib gspread
import gspread # https://docs.gspread.org/en/v6.0.0/user-guide.html
from google.oauth2.service_account import Credentials
```

```
scopes = ['https://www.googleapis.com/auth/spreadsheets']
credentials = Credentials.from_service_account_file('credentials.json', scopes=scopes)
client = gspread.authorize(credentials) # Use gspread to open the Google Sheet
spreadsheet = client.open_by_key("15R3QuyUFN_oVztI_P4lM9C5mjsYC7ZKUO-BogCdt81s") # Open the Google Sheets document by name
list_of_ampl_variables = [item[0] for item in ampl.get_variables()] # Retrieve a list of AMPL variables from the AMPL environment
variables_to_exclude = ['X', 'Cracking_X', 'Distill_X', 'Lubricating_X', 'Reform_X'] # List of variables to exclude from the export
filtered_ampl_variables = [var for var in list_of_ampl_variables if var not in variables_to_exclude]
for key_ampl in filtered_ampl_variables: # Iterate through each AMPL variable
ampl_dict = ampl.getVariable(key_ampl).getValues().toDict() # Retrieve AMPL variable values and convert them to a dictionary
filtered_dict = { # Filter the dictionary to exclude certain variables and include only values greater than 0.0001
key: value
for key, value in ampl_dict.items()
if value > 0.0001}
#if filtered_dict: # Only process further if there's data to be exported
if key_ampl not in [sheet.title for sheet in spreadsheet.worksheets()]: # Check if the current variable's name is not already a sheet in the spreadsheet
new_sheet_name = f"{key_ampl}" # Create a new sheet with a unique name based on the variable
spreadsheet.add_worksheet(title=new_sheet_name, rows=100, cols=20) # Add a new worksheet to the spreadsheet with a default size of 100 rows and 20 columns
sheet = spreadsheet.worksheet(key_ampl) # Get the worksheet object corresponding to the current variable
sheet.clear() # Clear the existing data in the sheet
first_key = next(iter(filtered_dict.keys()), None) # Get the first key to determine the data structure
if isinstance(first_key, tuple): # Check if the first key is a tuple, indicating multi-dimensional data
number_of_columns = len(first_key)
indices = [f'index{i}' for i in range(number_of_columns)]
data = [] # Prepare data for DataFrame
for key, value in filtered_dict.items():
data.append([str(k) for k in key] + [value]) # Convert each part of the key tuple to a string and append the value
columns = indices + [key_ampl+".val"] # Define column names including value column
else:
# Handle one-dimensional data
indices = [''] # Single index column
data = [[str(key), value] for key, value in filtered_dict.items()] # Convert key to string and pair it with the value
columns = indices + [key_ampl+".val"] # Define column names including value column
df = pd.DataFrame(data, columns=columns) # Create a DataFrame from the prepared data
data_list = [df.columns.tolist()] + df.values.tolist() # Convert DataFrame to list of lists
sheet.update(data_list,'A1') # Update the Google Sheet with the prepared data starting from cell A1
```

## Load Data into the Power BI Report from Google Sheets¶

To ensure we have the most up-to-date information, we’ll refresh the data in the Oil Refinery report. Once the data is updated, the latest model solution will be visualized on the “Result” page, providing a clear and accurate reflection of the current outcomes.