%%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{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_Cap_Max{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_Cap_Max{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, POLLUT} ; # Pollutant emissions
param crack_Cost{CR_MODE, DIST} >= 0 ; # Cost of process
param crack_Waste_Cost{CR_MODE, DIST_CR} >= 0 ; # Residue disposal cost
param crack_Cap_Max{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_Cap_Max{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} >= 0 ; # Step values for price elasticity
param demand_nStep_Value{1..nStep+1} >= 0 ; # 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]) ;
## 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);
## 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] * price_nStep_Value[n] - 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)^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[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_Cap_Max [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_Cap_Max[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_Cap_Max[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_Cap_Max[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 ;
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}:
Demand[p,s,t,n] >= demand_nStep_Value[n] * X[p,s,t,n] * seasonal_Base_Demand[p,t];
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] * demand_nStep_Value[n+1] * seasonal_Base_Demand[p,t];
## 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;