# Presentation: 122323

## Short-Term Scheduling under Uncertainty

- Marianthi Ierapetritou
- Department Chemical and Biochemical Engineering
- Piscataway, NJ 08854-8058

## Objective
Identify and  reduce bottlenecks at different levels 
Integration of the whole decision-making process

**Identify and  reduce bottlenecks at different levels **
**Integration of the whole decision-making process**

- Online Control

- Opportunity for Optimization

**Process Operations Decision Making**

**Notes:**
When we look at D-M Pro, the decision closest to the operation is online control which needs intermediate information acquisition and reaction. The next level is sh-t sch, which provides schedule from hours to days to the plant. When we move to weeks to months plan, we are looking at production planning problem. On top of the tree, it is supply chain management problem which considers multiple sites and multiple periods.
Time horizon for these problems are backward while the complexity and uncertainty increase along the structure arrows. The opportunity for optimization also grows because intuitive and experience couldn’t make right decisions when one is dealing with hundreds of processes and thousands of products. 
Our objective for this proposal is to identify and get over the bottlenecks at different levels, therefore integrate the whole D-M processes.  

## Short-term scheduling
Uncertainty in product prices, product demands, raw material availability, machine availability, processing times

Production Planning		
Longer time horizon under consideration (several months)
Larger number of materials and products
Uncertainty in facility availability, product demands, orders, raw materials

Supply chain management
Multiple sites involving production, inventory management, transportation 
Longer planning time horizon (couple of years)
Uncertainty in material availability, costs, transportation

- Uncertainty in product prices, product demands, raw material availability, machine availability, processing times
- Production Planning
- Longer time horizon under consideration (several months)
- Larger number of materials and products
- Uncertainty in facility availability, product demands, orders, raw materials
- Supply chain management
- Multiple sites involving production, inventory management, transportation
- Longer planning time horizon (couple of years)
- Uncertainty in material availability, costs, transportation

**Uncertain Parameters**

**Notes:**
As mentioned before, these challenges for short-term scheduling problem are uncertainties in pro....  Lar because of dozens of ...... which are translated to an computational expensive model. I will get into detail in later slides.
For production planning, challenges in large scale also exist because even longer time horizon is considered here as well as larger number of materials and processes. Uncertainty also arises in facility availability, which product in demand in future and so on.
Supply chain management involves production planning, inventory management, transportation management and other processes that needs to be optimized. How to integrate them together and contribute to the overall profit is the goal of supply chain optimization as we will talk about at last.

## Process Plant				Optimal Schedule

Given:					       	Determine:
Raw Materials, Required Products, 	       	Task Sequence, 
Production Recipe, Unit Capacity	       	Exact Amounts of material
							Processed
Scheduling objectives :		       		
Economic 	Maximize Profit, Minimize Operating Costs, 
			Minimize Inventory Costs
Time Based 	Minimize Makespan, Minimize Tardiness

- Given:					       	Determine:
- Raw Materials, Required Products, 	       	Task Sequence,
- Production Recipe, Unit Capacity	       	Exact Amounts of material
- 							Processed
- Scheduling objectives :
- Economic 	Maximize Profit, Minimize Operating Costs,
- 			Minimize Inventory Costs
- Time Based 	Minimize Makespan, Minimize Tardiness

**Short-term Scheduling**

## Binary variables to allocate tasks to resources

Continuous variables to represent  timing and material variables 

Mixed Integer Linear Programming Models

Smaller models that are  computationally efficient and tractable

- Continuous variables to represent  timing and material variables
**M**ixed **I**nteger **L**inear **P**rogramming Models
- Smaller models that are  computationally efficient and tractable

**Continuous Time Formulation**

## Deterministic Scheduling Formulation

- minimize H   or  maximize price(s)d(s,n)
- subject to 	wv(i,j,n)  1
- st(s,n) = st(s,n-1) – d(s,n) + **P**b(i,j,n-1) + **c**b(i,j,n)
- st(s,n)  stmax(s)
- Vmin(i,j)wv(i,j,n)  b(i,j,n)  Vmax(i,j)wv(i,j,n)
- d(s,n)  r(s)
- Tf(i,j,n) = Ts(i,j,n) + (i,j)wv(i,j,n) + (i,j)b(i,j,n)
- Ts(i,j,n+1)  Tf(i,j,n) – U(1-wv(i,j,n))
- Ts(i,j,n)  Tf(i’,j,n) – U(1-wv(i’,j,n))
- Ts(i,j,n)  Tf(i’,j’,n) – U(1-wv(i’,j’,n))
- Ts(i,j,n)  H, Tf(i,j,n)  H

- Duration Constraints

- Demand Constraints

- Allocation Constraints

- Capacity Constraints

- Material Balances

- Objective Function

- n

- s

**(i,j)**

**M.G.Ierapetritou and C.A.Floudas. *****Effective continuous-time formulation for short-term scheduling. 1. Multipurpose batch processes*****. 1998**

**Notes:**
Before going to the details of these two approaches, I’ll first introduce the basic model for the case of deterministic scheduling. The objective is normally to maximize the profit or minimize the makespan. The constraints mainly include allocation constraints, material balance constraints, capacity, demand and duration constrains. Binary variable wv denotes the assignment of task i to unit j at event point n. So this model is a mixed integer linear programming problem.

## Increased Complexity: Parameter Fluctuations

**Increased Complexity: Parameter Fluctuations**

**Notes:**
This graph shows the annual price change of crude oil. It will be interesting to investigate effects of oil price fluctuation, and other uncertainties on refinery operations. 

## Uncertainty in Short-Term Scheduling

- Price of P1 is an uncertain parameter. Considering time horizon of 16 hours, $1 increase results in the following different production schedules.

**Uncertainty impacts the optimal schedule**

**Uncertainty in Short-Term Scheduling**

## Uncertainty in Short-Term Scheduling

**0**

**8**

**2**

**4**

**6**

**1**

**3**

**5**

**7**

**0**

**8**

**2**

**4**

**6**

**1**

**3**

**5**

**7**

**50.00**

**50.00**

**10.40**

**64.60**

**10.03**

**50.93**

**64.04**

**4.63**

**50.93**

**55.56**

**74.07**

**4.63**

**50.93**

**44.44**

**74.07**

**50.93**

**55.56**

- Deterministic Schedule

- Robust Schedule

- heating

- heating

- reaction 1

- reaction 2

- reaction 3

- reaction 2

- reaction 1

- reaction 2

- reaction 3

- separation

- reaction 2

- reaction 2

- reaction 3

- reaction 3

- reaction 1

- reaction 1

- separation

- demand (product 2) = 50*(1 _+_ 60%)

- demand (product 2) = 50

- Standard Deviation = _0.29_

- E(makespan) = _7.24hr_

- E(makespan) = _8.15hr_

- Standard Deviation = _2.63_

**Uncertainty in Short-Term Scheduling**

**Notes:**
By solving this optimization problem, we’ll get an optimal schedule under nominal data. For this example, the nominal demand is 50 and this is the deterministic schedule with makespan equal to 7 hours. However, in reality, the demand could have a deviation of plus minus 60 percent. In this case, the lower schedule can accommodate the demand fluctuation much better than the deterministic schedule even it has a longer makespan when demand is 50 but it . 

## Literature Review: Representative Publications

-  Reactive Scheduling

-  S.J.Honkomp, L.Mockus, and G.V.Reklaitis.  *A framework for schedule evaluation with processing uncertainty*. Comput. Chem. Eng. 1999, 23, 595
-  J.P.Vin and M.G.Ierapetritou. *A new approach for efficient rescheduling of multiproduct batch plants*. Ind. Eng. Chem. Res., 2000, 39, 4228

- Handles uncertainty by adjusting a schedule upon realization of the uncertain parameters or occurrence of unexpected events

-  Stochastic Programming

- Uncertainty is modeled through discrete or continuous probability functions

-  J.R.Birge and M.A.H.Dempster.  *Stochastic programming approaches to stochastic scheduling*. J. Global. Optim. 1996, 9, 417
-  J.Balasubramanian and I.E.Grossmann. *A novel branch and bound algorithm for scheduling flowshop plants with uncertain processing times*. Comput. Chem. Eng. 2002, 26, 41

## Literature Review: Representative Publications

**Literature Review: Representative Publications**

-  Robust Optimization

-  X.Lin, S.L.Janak, and C.A.Floudas.  *A new robust optimization approach for scheduling under uncertainty – I. bounded uncertainty*. Comput. Chem. Eng. 2004, 28, 2109

- Produces “robust” solutions that are immune against uncertainties

-  Fuzzy Programming

-  H.Ishibuchi, N.Yamamoto, T.Murata and Tanaka H. *Genetic algorithms and neighborhood search algorithms for fuzzy flowshop scheduling problems *. Fuzzy Sets Syst. 1994, 67, 81
-  J.Balasubramanian and I.E.Grossmann. *Scheduling optimization under uncertainty- an alternative approach*. Comput. Chem. Eng. 2003, 27, 469

- Considers random parameters as fuzzy numbers and the constraints are treated as fuzzy sets

-  MILP Sensitivity Analysis

-  Z.Jia and M.G.Ierapetritou. *Short-term Scheduling under Uncertainty Using MILP Sensitivity Analysis.*  Ind. Eng. Chem. Res. 2004, 43, 3782

- Utilizes MILP sensitivity analysis methods to investigate the effects of
- uncertain parameters and provide a set of alternative schedules

## Uncertainty in Scheduling

**Uncertainty in Scheduling**

## Preventive Scheduling

- model robustness

- solution robustness

- Data perturbation

-  New alternative schedules

- Deterministic schedule

- LB/UB on objective function

- MILP sensitivity analysis framework

- Robust optimization
- method

- A set of solutions represent trade-off between various objectives

**Preventive Scheduling**

**Notes:**
The formulations of all these scheduling problems are based on continuous time representation, which means an event point is used only if a task takes place. However, in discrete time formulation, the time horizon is divided into equal sized intervals beforehand, so extra points are used if nothing happens. The main advantage of continuous time representation is the model size can be much smaller, which will significantly reduce the computational time. For scheduling under uncertainty, we proposed two approaches. One is a framework based on MILP sensitivity analysis and the other is to develop a robust optimization model which considers  both model robustness and solution robustness. A solution is...

## Preventive Scheduling

_**Robust Optimization**_

_**MILP Sensitivity Analysis**_

- minimize H   or  maximize price(s)d(s,n)
- subject to 	wv(i,j,n)  1
- st(s,n) = st(s,n-1) – d(s,n) + **P**b(i,j,n-1) + **c**b(i,j,n)
- st(s,n)  stmax(s)
- Vmin(i,j)wv(i,j,n)  b(i,j,n)  Vmax(i,j)wv(i,j,n)
- d(s,n)  r(s)
- Tf(i,j,n) = Ts(i,j,n) + (i,j)wv(i,j,n) + (i,j)b(i,j,n)
- Ts(i,j,n+1)  Tf(i,j,n) – U(1-wv(i,j,n))
- Ts(i,j,n)  Tf(i’,j,n) – U(1-wv(i’,j,n))
- Ts(i,j,n)  Tf(i’,j’,n) – U(1-wv(i’,j’,n))
- Ts(i,j,n)  H, Tf(i,j,n)  H

- Mixed-integer Linear Programming

**Notes:**
Since the scheduling problems are modeled as mixed integer linear programming problems, our first approach is based on MILP sensitivity analysis.

## Questions to Address

- mixing

- reaction

- separation

- H (time horizon)

**10**

**0**

**2**

**4**

**8**

**55**

**55**

**50**

- mixing

- reaction

**15**

**6**

**15**

**20**

- separation

- Can the schedule accommodate the demand fluctuation?

- How the capacity of the units affect the production objective?

** **What is the effect of processing time at the objective value?

**Questions to Address**

**Notes:**
Through this approach, we would like to address questions such as what is the effect of processing time at the objective, how the capacity of the units affects the production, and can the schedule handle the demand fluctuation. 

## Inference-based MILP Sensitivity Analysis

**Inference-based MILP Sensitivity Analysis**

- **i****P **A**ij**_u_**j****P** + s**j**(u**j** – _u_**j**) - **i**a**i**  r**P**
- s**j****P**  **i****P**A**ij**, s**j****P**  -q**j****P**, j = 1,...,n
- r**P** = -q**j****P**u**j****P** + **P**a – z**P** +z**P**

- c**j**_u_**j****P** - s**j****P**(u**j****P** – _u_**j****P**)  -r**P**
- s**j****P**  -c**j**, s**j****P**  -q**j****P**, j = 1,...,n
- q**j****P **= **i****P**A**ij** - **i****P**c**j**

- - for the perturbations A and a

- - for the perturbations c

- Bound z  z***** - z holds if there are s**1****P****,...,**s**n****P**** **that satisfy:

- minimize   z = cx
- subject to   Ax  a
- 0 x  h, xj integer, j=1,...k

- minimize   z = (c + c)x
- subject to   (A + A)x  a + a
- 0 x  h, xj integer, j=1,...k

- Aim: Determine under what condition z  z***** - z remains valid

- _Partial assignment at node p_

- *M.W.Dawande and J.N.Hooker, 2000

**Notes:**
In their work, it is found that if a minimization MILP problem is perturbed, the bound that the objective decreases at most delta z is valid as long as these two linear systems are satisfied at each leaf node in the b&b tree. Leaf nodes are nodes at which there’s no further requirement of branching. This term denotes the partial assignment of variable x imposed on node p by branching cuts. Comparing with most existing approaches, this analysis is more general, in that it allows perturbations of any data, including the objective function and constraints coefficients, as well as the right hand sides.

## Proposed Uncertainty Analysis Approach

- Solve the deterministic scheduling problem using B&B tree

**Proposed Uncertainty Analysis Approach**

-  Range of objective change for certain parameter change

- MILP Sensitivity Analysis

## Robustness Estimation

- Obtain sequence of tasks
- from original schedule

- Generate random demands in expected range

- Makespan minimization is considered as the objective

- Makespan to meet a particular demand is found using the sequence of tasks derived from original schedule
- Binary variables corresponding to allocation of tasks are fixed
- Batch sizes and Starting and Finishing times of tasks are allowed to vary

**Robustness Estimation**

## Robustness under Infeasibility

- Corrected Standard Deviation:

-  H*act*  =  H*p *      if scenario is feasible
-          =  H*corr*   if the scenario is infeasible

**Robustness under Infeasibility**

**J.P.Vin and M.G.Ierapetritou. *****Robust short-term scheduling of multiproduct batch plants under demand uncertainty*****. 2001**

## Case Study 1

**Case Study 1**

- S1

- S2

- S3

- S4

- wv(i1,j1,n0)

- wv(i1,j1,n1)

- wv(i2,j2,n1)

- wv(i3,j3,n2)

- wv(i2,j2,n2)

**3.0**

**5.17**

**3.0**

**7.65**

**5.83**

**7.16**

**infeasible**

**8.14**

**10.16**

**7.16**

**9.83**

**infeasible**

**infeasible**

**8.33**

**9.87**

**8.83**

**7.16**

**8.83**

**9.98**

**infeasible**

**9.83**

**(Schedule 1)**

**0**

**0**

**0**

**0**

**0**

**0**

**0**

**0**

**0**

**0**

**1**

**1**

**1**

**1**

**1**

**1**

**1**

**1**

**1**

**1**

- _Effect of demand _ d~[20, 100]

- -0.097 d  H

- dnom = 50

- H**’**  Hnom + 0.097d = 12.73h

- d**’** = 80

- Hnom = 9.83h

- mixing

- reaction

- purification

- B&B tree with nominal demand

**Notes:**
Our case study involves a production line consisting of a mixer, a reactor and a purifier. The demand is considered to be uncertain varying within the range of 20 to 100. The problem is first solved at nominal demand 50 and the optimal schedule is schedule 1 with makespan 9.83 hours. The value next to each node is the objective value of the relaxed LP and these are the binary variables being branched at each level Applying the sensitivity analysis, it is found that if the demand is increased by delta d, the new makespan is increased by at most 0.097 delta d. When the demand is increased to 80, the current schedule becomes infeasible. Using this inequality, the new makespan is less than 12.73 hours.

## Case Study 1

**Case Study 1**

**(12.13)**

- schedule 1

- schedule 2

- schedule 3

- H**nom**(h)

- H**avg**(h)

- SD**corr**

- 9.83

- 14.20

- 5.52

- 10.77

- 11.56

- 1.61

- 10.91

- 11.79

- 2.17

- wv(i1,j1,n0)

- wv(i1,j1,n1)

- wv(i2,j2,n1)

- wv(i3,j3,n2)

- wv(i3,j3,n3)

- wv(i2,j2,n2)

**3.0**

**5.17**

**3.0**

**7.65**

**5.83**

**7.16**

**infeasible**

**8.14**

**10.16**

**7.16**

**9.83**

**infeasible**

**infeasible**

**8.33**

**9.87**

**8.83**

**7.16**

**8.83**

**9.98**

**infeasible**

**9.83**

**(Schedule 1)**

**0**

**0**

**0**

**0**

**0**

**0**

**0**

**0**

**0**

**0**

**1**

**1**

**1**

**1**

**1**

**1**

**1**

**1**

**1**

**1**

**(12.73)**

**(17.97)**

****

****

****

- Schedule Evaluation

**Notes:**
So we solve the relaxed LP of the leaf nodes with demand equal to 80, and only these two nodes have objective less than or equal to 12.73. We continue branching on these two nodes and find optimal schedule 2 and feasible schedule 3. The schedules are then evaluated based on nominal makespan, average makespan, and corrected standard deviation. Although schedule 1 provides the best solution at nominal point, it has poor performance over the entire interval. Schedule 2 has the highest robustness and shortest average makespan, and schedule 3 provides a compromise solution with robustness between the ones of schedule 1 & 2, and nominal average makespan close the schedule 2. 

## Case Study 1

**schedule 1 (optimal when d ****≤ 50)**

**schedule 2 (optimal when d ****≥**** 50)**

**schedule 3**

**Case Study 1**

**Notes:**
This is the graphic presentation of the results. We can see that although schedule 1 has the shortest makespan when demand is less than 50, schedule 2 and 3 have more flexibility in the entire demand range.

## Case Study 1

- _Effect of processing time_  T(i1,j1) ~ [2.0, 4.0]

- profit**’**** ** profitnom + 24.48T = 47.04

- Tnom = 3.0

- T**’** = 4.0

- profitnom = 71.52

- schedule 1

- schedule 2

- schedule 3

- profit**nom**

- profit**avg**

- SD**corr**

- 71.52

- 66.98

- 26.9

- 65.27

- 9.33

- 65.27

- 65.17

- 10.49

- wv(i1,j1,n0)

- wv(i1,j1,n1)

- wv(i2,j2,n1)

- wv(i3,j3,n2)

- wv(i3,j3,n3)

- wv(i2,j2,n2)

**100**

**100**

**50**

**100**

**100**

**96.05**

**78.42**

**50**

**(75)**

**0**

**0**

**0**

**0**

**0**

**0**

**0**

**1**

**1**

**1**

**1**

**1**

**1**

**1**

**71.52**

**100**

**50**

**75**

**(Schedule 1)**

**50**

**75**

**(75)**

**1**

**0**

**50**

**72.46**

**(62.11)**

**78.42**

- 64.61

**Notes:**
There are other uncertainties in scheduling problems, processing time is one of them. With the similar procedure, we can discover 3 schedules when the processing time of the mixer varies from 3 to 4 hours.

## Shortcoming of Proposed Approach

- Since the entire analysis is based on a single tree among a large number of possible branch-and-bound trees that can be used to solve the MILP, it provides conservative sensitivity ranges.

## Parametric Programming

**Parametric Programming**

- z() = min  cTx + dTy
- subject to   Ax + Dy  b
- xL  x  xU
- L    U
- x  Rm, y  (0,1)T

- z() = min  cTx + dTy
- subject to   Ax + Dy - r b
- cTx + dTy – z0 -  = 0
- yi - yi  F1 - 1
- xL  x  xU
- L’    U’
- x  Rm, y  (0,1)T

- iF1

- iF 0

- b = b0 + r

- LP sensitivity analysis:

- z() = z**0 **+ 
- L  L’    U’  U

- solved at b = b0+**L**r
- optimal solution (x*****,y*****)

- Integer cut to exclude current optimal solution

-  b****[b0+**L**r, b0+**U**r]

- break point ’, new
-   optimal solution (x*****’, y*******’**)

- Fix integer variables at y*****

- A.Pertsinidis et al. Parametric optimization of MILP programs and a framework for the parametric optimization of MINLPs. 1998

**Notes:**
Parametric programming looks for all the optimal schedule in a given parameters range. In this MILP problem, right hand side vector b varies with respect of theta. The problem is first solved at one point and optimal solution is found. Then LP sensitivity analysis is performed with fixed integer values. The result and an integer cut is added to the original problem. By solving this new MILP problem, we can find the break point theta prime, and the successor optimal solution. Repeating this procedure, we can get all optimal solution in the range of theta lower to theta upper.

## Multiparametric MILP

-  J. Acevedo and E.N.Pistikopoulos. *A Multiparametric Programming Approach for Linear Process Engineering Problems under Uncertainty*. 1997

- Select a branching variable

- Solve the fully relaxed problem

- Solve the mpLP at the nodes

- Compare the solution with the current UB, update the optimal function in the uncertain space

- Based on simplex algorithm, check the neighboring bases of the LP tableau

## Shortcomings of Existing Approach

-  mpLP approach requires retrieving the LP tableaus and visiting the neighbor bases

-  Solve mpLP at every node in the B&B tree during the branch and bound procedure: can be a computationally expensive effort

## Proposed Analysis on the RHS for MILPs

- minimize   z = cx
- subject to   Ax  a
-  x  0, xj Є (0, 1), j=1,...k

- minimize     z = cx
- subject to   Ax  a + a
- x  0, xj Є (0, 1) , j=1,...k

- Develop a framework to investigate the effect of Δa on the optimal solution x and objective value z

-  A set of optimal integer solutions
-  Critical regions
-  Optimal functions

## Proposed Approach: Single Uncertain Parameter

- Find Δamax that leaves the structure of the tree unchanged

- Solve the original problem at the nominal value using a branch and bound method

- Collect zp, λp at each leaf node p

- For a = Δamax + ε
- update the B&B tree

## Determine Δamax and Update the B&B Tree

- Case 1: 	A descent node of the node 0

- Case 2: A descent node of other leaf node

- Case 3: Node 0, but the basis has changed

- _The new optimal node can be:_

**Determine ****Δ****a****max**** and Update the B&B Tree**

- P

- Update the B&B tree at a = Δamax +ε

- where node 0  is the optimal node

**0**

## Proposed Approach: Multiple Uncertain Parameters

- mpLP at the leaf nodes

- mpLP algorithm

-     Compare the critical regions
-     with the current upper bounds
-     Update the B&B tree

- Solve the original problem at the nominal value using a branch and bound method

## mpLP Algorithm at the Leaf Nodes

- maximize	cx* – z
- subject to   z = max{z(k) + λ(k)θa + β(k)θb}
-  	          a0 ≤ θa ≤ a0 + Δa
- 	          b0 ≤ θb ≤ b0 + Δb

- Bilevel  Linear Programming

- At each iteration, solve

- where x* = argmin cx
- subject to  Ax ≥ θ

- θ

- cx*

- max{z(k) + λ(k)θ}

- cx*

- current optimal functions

## mpLP Algorithm at the Leaf Nodes

- maximize   {min cx| Ax  θ} – z
- subject to   z  z(k) + λ(k)θa + β(k)θb
-  	          a0 ≤ θa ≤ a0 + Δa
- 	          b0 ≤ θb ≤ b0 + Δb

- current optimal functions

- check if there is any point at which max{z(k) + λ(k)θa + β(k)θb} is less than cx*

- Stop when the objective = 0

- Include an additional constraint z  z(k+1) + λ(k+1)θa + β(k+1)θb  to above problem

## Solve the Bilevel Programming Problem

- max  	    θy
- s.t.   	ATy ≤ c
-         a0 ≤ θa ≤ a0 + Δa
-         b0 ≤ θb ≤ b0 + Δb
-    	Y  0

- min   	cx
- s.t.     Ax  θ
-        a0 ≤ θa ≤ a0 + Δa
-        b0 ≤ θb ≤ b0 + Δb
- 	 x  0

- Replace the inside optimization problem by its dual

- max  	    θy - z
- s.t.   	ATy ≤ c
- z  z(k) + λ(k)θa + β(k)θb
-         a0 ≤ θa ≤ a0 + Δa
-         b0 ≤ θb ≤ b0 + Δb
-    	Y  0

- Single optimization problem

- Convert the relaxed LP to its dual form

- Strong Duality Theorem

- Bilinear objective
- Linear constraints

- Solve with global optimization solver BARON

## Compare the Optimal Functions of the Leaf Nodes

- z1UB = z1*UB + λ1UBθa + β1UBθb

- 1

- 2

- z2(2) = z2*(2) + λ2(2)θa + β2(2)θb

- CR2(2)

- CR2(1)

- CR2(3)

- CR1(1)

- CR1UB

- CR1UB   CR2(2) = CRint

## Compare the Optimal Functions of the Leaf Nodes

- min 	    Є
- s.t.   	z1UB = z2(2) + Є
- z1UB = z1*UB + λ1UBθa + β1UBθb
- z2(2) = z2*(2) + λ2(2)θa + β2(2)θb
- 	θa,θb Є CRint

- Redundancy test on constraint z1UB ≥ z2(2)

- _Case 1_: 	Problem is infeasible: z1UB is smaller in CRint

- _Case 2_: Є > 0: the constraint is redundant. z2(2) is smaller in CRint  The optimal function is updated to be z2(2) if node 2 is an integer node, otherwise, do not update.

- Compare optimal function z1UB and z2(2) in region CRint

- _Case 3_: Є < 0: the constraint is not redundant. CRint is divided into two parts. z1UB is smaller on one side and z2(2) is smaller on the other side. The two regions are divided by z1UB ≤ z2(2).

## Advantage of Proposed Approach

-  The new mpLP approach can efficiently determine the optimal functions and critical regions without retrieving the LP tableaus and visiting the neighbor bases

-  Solve mpLP at only the leaf nodes in the B&B tree instead of every node during the branch and bound procedure – reduce the computational efforts significantly

## Case Study: Single Uncertain Parameter 

- min    z = 2x1 + 3x2 + 1.5x3 + 2x4 + 0.5x5
- s.t.   	2x1 + x2 + x3 ≥** **7 + Δa
- 	2x2 + x4 + x5 ≥ 4
- 	x3 + x4 – x5 ≥** **0
- 	2x1 – x2 – x3 + x5 ≥ 4
-         1≤ x ≤ 3, xj Є (0,1), j = 3,4,5

- _Step 1:_

**3**

**2**

** ****x****3**** = 1**

** ****x****3**** = 0**

** ****x****4**** = 1**

** ****x****4**** = 0**

**0**

**12**

**11.5***

**12**

**(0,1,1)**

- _Step 2:_

- Linear sensitivity analysis on node 0 -- Δamax = 0

- _Step 3:_

- For Δa = Δamax + ε, node 0 yields noninteger solution. Update the B&B tree.

-  Δa = 0

## Case Study: Single Uncertain Parameter 

**3**

**2**

** ****x****3**** = 1**

** ****x****3**** = 0**

** ****x****4**** = 1**

** ****x****4**** = 0**

**0**

**12**

**11.8***

**12**

**λ**** = 3**

** ****x****5**** = 1**

** ****x****5**** = 0**

**12.1**

**3**

**2**

** ****x****3**** = 1**

** ****x****3**** = 0**

**1**

**0**

**0**

**12.2**

**12.1**

**12***

**1**

**0**

**12.2**

**1**

**1**

**0**

**0**

**12***

**4**

**5**

**λ**** = 0**

**λ**** = 0**

**λ**** = 1**

**λ**** = 1**

**λ**** = 2**

- _Step 2:_

- Δamax = min{Δabasis, min{             } }

**P**

- = min{ 2, min {                       }}

- =

- _Step 3:_

- For Δa = Δamax + ε, node 0 is intersected by node 2 & 3. Update the B&B tree

- _Step 3:_

- _Step 2:_

**...**

**...**

**Case Study: ****Single Uncertain Parameter**** **

## Case Study: Single Uncertain Parameter 

- 0

- 0.2

- 0.4

- 0.6

- 0.8

- 1.0

- 1.2

- 1.4

- 1.6

- 1.8

- 2

- 11.5

- 12

- 12.5

- 13

- 13.5

- 14

**(0,1,1)**

**(1,0,1) & (0,0,0)**

**(1,0,1)**

**(1,0,1)**

- Final optimal solution:

- z*

- Δa

**Case Study: ****Single Uncertain Parameter**** **

## Case Study: Multiple Uncertain Parameters 

- min    z = -3x1 - 8x2 + 4y1 + 2y2
- s.t.   	x1 + x2 ≤** **13 + θ1
- 	5x1 - 4x2 ≤ 20
-  	-8x1 + 22x2 ≤ 121 + θ2
- 	4x1 + x2 ≥** **8
- 	x1 – 10y1 ≤ 0
- 	x2 – 15y3 ≤ 0
- x ≥ 0, y Є (0,1), 0 ≤ θ1,θ2 ≤ 10

- _Step 1:_

**1**

** ****y****1**** = 0**

** ****y****1**** = 1**

** ****y****2**** = 0**

** ****y****2**** = 1**

**2**

**-8**

**-70.5***

- (θ1, θ2) = (0, 0)

## Case Study: Multiple Uncertain Parameters 

- s.t.   	-d1 – 5d2 + 8d3 + 4d4 – d5 ≤** **-3
- 	-d1 + 4d2 – 22d3 + d4 – d6 ≤ -8
- 		10d5 – d7 ≤** **4
- 		15d6 – d8 ≤ 2
- 	z ≥ -70.5 – 4.3333θ1 – 0.1667θ2
-        		0 ≤ θ1, θ2 ≤ 10

- _Step 2:_

- mpLP on the leaf nodes

- Node 1

- z1(1) = -70.5 – 4.3333θ1 – 0.1667θ2     at (θ1, θ2) = (0, 0)

- obj = 16.74    (θ1, θ2) = (10, 0)

- z1(2) = -97.0909 – 0.1667θ2     at (θ1, θ2) = (10, 0)

**Case Study: ****Multiple Uncertain Parameters**** **

## Case Study: Multiple Uncertain Parameters 

- s.t.   	-d1 – 5d2 + 8d3 + 4d4 – d5 ≤** **-3
- 	-d1 + 4d2 – 22d3 + d4 – d6 ≤ -8
- 		10d5 – d7 ≤** **4
- 		15d6 – d8 ≤ 2
- 	z ≥ -70.5 – 4.3333θ1 – 0.1667θ2
- 	z ≥ -97.0909 – 0.3636θ2
-        		0 ≤ θ1, θ2 ≤ 10

- obj = 0, stop

- z1(1) ≥ z1(2)            0.07333θ1 – 0.00333θ2 ≤ 0.45

- The two critical regions are divided by

**Case Study: ****Multiple Uncertain Parameters**** **

## Case Study: Multiple Uncertain Parameters 

- Node 2

- z2 = 8

- CR2 =

- θ1 ≤ 10

- θ2 ≤ 10

- _Step 3:_

- Compare critical regions and determine optimal functions

- 0.07333θ1 – 0.00333θ2 ≤ 0.45

- θ2 ≤ 10

- z1(1) = -70.5 – 4.3333θ1 – 0.1667θ2

- CR1(1) =

- z1(2) = -97.0909 – 0.3636θ2

- 0.07333θ1 – 0.00333θ2 ≥ 0.45

- θ1, θ2 ≤ 10

- CR1(2) =

- CR1(1) and CR2

- CR1(1) n CR2 = CR1(1)

**Case Study: ****Multiple Uncertain Parameters**** **

## Case Study: Multiple Uncertain Parameters 

- s.t.   	-97.0909 – 0.3636θ2 + Є =** **-8
- 	 0.07333θ1 – 0.00333θ2 ≥ 0.45
-        		0 ≤ θ1, θ2 ≤ 10

- Є > 0,

- z1(2) ≤ z2  is redundant in CR1(2)

- min   		    Є
- s.t.   	-70.5 – 4.3333θ1 - 0.1667θ2 + Є =** **-8
- 	 0.07333θ1 – 0.00333θ2 ≤ 0.45
-        		 θ2 ≤ 10

- Є < 0,

- z1(1) ≤ z2

- CR1(2) and CR2

- CR1(2) n CR2 = CR1(2)

- CR1(1)

- CR1(2)

**Case Study: ****Multiple Uncertain Parameters**** **

## Case Study: Multiple Uncertain Parameters 

- 0.07333θ1 – 0.00333θ2 ≤ 0.45

- θ2 ≤ 10

- z1(θ) = -70.5 – 4.3333θ1 – 0.1667θ2

- CR1(1) =

- z2(θ) = -97.0909 – 0.3636θ2

- 0.07333θ1 – 0.00333θ2 ≥ 0.45

- θ1, θ2 ≤ 10

- CR1(2) =

- y* = (1, 1)

- Final optimal solution:

**Case Study: ****Multiple Uncertain Parameters**** **

## Uncertainty Analysis on the Objective Function Coefficients 

- minimize   z = cx
- subject to   Ax  θ
-  x  0, xj Є (0, 1), j=1,...k

- minimize     z = (c + Δc)x
- subject to   Ax  θ
- x  0, xj Є (0, 1) , j=1,...k

- z* = max{z(k) + λ(k)θa + β(k)θb}

- Unlike the case of uncertain RHS, where the optimal objective value

- Here,  z* = min{z(k) + λ(k)c1 + β(k)c2}

## mpLP Algorithm at the Leaf Nodes

- subject to   	Ax ≥ θ
- 	         z ≤ z(k) + λ(k)c1 + β(k)c2
-  	          c10 ≤ c1 ≤ c10 + Δc1
- 	          c20 ≤ c2 ≤ c20 + Δc2

**mpLP Algorithm at the Leaf Nodes**

- maximize	z – (c + Δc)x*
- subject to   z = min{z(k) + λ(k)θa + β(k)θb}
-  	          c10 ≤ c1 ≤ c10 + Δc1
- 	          c20 ≤ c2 ≤ c20 + Δc2

- Bilevel  Linear Programming

- where x* = argmin (c+Δc)x
- subject to  Ax ≥ θ

- One level  NLP

## Uncertainty Analysis on the Constraint Coefficients

- minimize   z = cx
- subject to   Ax  θ
-  x  0, xj Є (0, 1), j=1,...k

- minimize     z = cx
- subject to   (A + ΔA)x  θ
- x  0, xj Є (0, 1) , j=1,...k

- maximize	cx* – z
- subject to   z = max{z(k) + λ(k)a1 + β(k)a2}
-  	          a10 ≤ a1 ≤ a10 + Δa1
- 	          a20 ≤ a2 ≤ a20 + Δa2

- Bilevel  Linear Programming

- where x* = argmin cx
- subject to  Ax ≥ θ

- At mpLP procedure, need to solve

## Solve Bilevel Linear Programming Problem

- The most popular method is “Kuhn-Tucker” approach

- min		F(x, y) = c1x + d1y
- subject to   	A1x + B1y ≤ b1

- min		f(x, y) = c2x + d2y
- subject to   	A2x + B2y ≤ b2

**x****Є****X**

**y****Є****Y**

- Replace with its KKT condition and add to the upper level problem

## Preventive Scheduling

_**Robust Optimization**_

_**MILP Sensitivity Analysis**_

- Expected Makespan/Profit

- Model Robustness

- Solution Robustness

- Objective =

**Notes:**
The approach I have presented so far does not consider robustness before solving for the alternative schedules. The schedules are found using the framework and then evaluated.  In our next approach, the robustness characteristics are included in the objective when we solve the scheduling optimization problems. 

## Robust Optimization

**Robust Optimization**

- _*__Upper Partial Mean_

-  S.Ahmed and N. Sahinidis. *Robust process planning under uncertainty*. 1998

- st**k**(s,n) = st**k**(s,n-1) – d**k**(s,n) + **P**(s,i)b**k**(i,j,n-1) + **c**b**k**(i,j,n)
- st**k**(s,n)  stmax(s)
- Vmin(i,j)wv(i,j,n)  b**k**(i,j,n)  Vmax(i,j)wv(i,j,n)
- d**k**(s,n) + **slack****k**(s)  r(s)
- Tf**k**(i,j,n) = Ts**k**(i,j,n) + (i,j)wv(i,j,n) + (i,j)b**k**(i,j,n)
- Ts**k**(i,j,n+1)  Tf**k**(i,j,n) – U(1-wv(i,j,n))
- Ts**k**(i,j,n)  H**k**, Tf**k**(i,j,n)  H**k**
******k ** H**k **–  P**k**H**k**, ******k ** 0

**k**

**n**

**(i,j)**

- minimize

- Average Makespan

- Model Robustness

- Solution Robustness

- _Unsatisfied Demand_

- subject to 	wv(i,j,n)  1

**Notes:**
This is our proposed robust optimization model. The superscript k denotes uncertain scenarios. The first objective is to minimize the average makspan over all scenarios. To consider model robustness, we include slack variables in the demand constraint and penalize the infeasible part in the objective function. To consider solution robustness, we use upper partial mean to measure objective variability. So instead of penalizing the all the deviations, we only penalize the positive deviations from the average makespan.

## Multiobjective Optimization

- _Pareto Optimal Solution:_

- A point x*єC is said to be Pareto optimal if and only if there is no such xєC that fi(x) ≤ fi(x*) for all i={1,2,...,n} , with at least one strict inequality.** **

- Min  F(x) =

- f1(x)

- f2(x)

- fn(x)

**:**

**:**

- xєC

- C = {x: h(x) = 0, g(x) ≤ 0, a ≤ x ≤ b}

**Notes:**
Now we need to solve a multiobjective optimization problem. Instead of solving for one optimal point, we are looking for Pareto optimal surface. Pareto optimal point basically means that from that point you cannot improve one objective without deteriorating the other one. For example, arc ACB is a Pareto surface.

## Normal Boundary Intersection (NBI)

**Normal Boundary Intersection (NBI)**

- f1*

- f2*

- F*

- Max  t

- g(x) ≤ 0

- h(x) = 0

- a ≤ x ≤ b

- x,t

- Φω+ t n = F(x) – F*

- ^

- s.t.

- (Utopia point)

- Min  F(x) =

- f1(x)

- f2(x)

- NBIω:

- A point in the Convex Hull of Individual Minima (CHIM)

-  I. Das and J. Dennis. *NBI: A new method for generating the Pareto surface in nonlinear multicriteria optimization problems*. 1996

- _Advantage:_ can produce a set of evenly distributed Pareto points independent of relative scales of the functions

**Notes:**
One approach of obtaining pareto surface is normal boundary intersection. The advantage of this technique is that it can produce a set of evenly distributed Pareto points independent of relatives scales of the functions. Assume we have two objectives, f1 and f2. First the individual minima points are determined by minimizing f1 or f2 only, which are point A and B. Then starting from the points on line segment AB, we search on the perpendicular direction towards the utopia point in the feasible region so as to find the Pareto surface, which is the solution of this subporblem. The .. represents a point on AB. 

## Case Study

- S1

- S2

- S3

- S4

- f1(x*) =

- 6.46
- 54
- 0.09

- f2(x*) =

- 11.5
- 0
- 0.66

- f3(x*) =

- 6.77
- 50
- 0

- mixing

- reaction

- purification

- F* =

- 6.46
- 0
- 0

- Φ =

- 0
- 54
- 0.09

- 4.14
- 0
- 0.66

- 0.31
- 50
- 0

- 0
- 54
- 0.09

- 4.14
- 0
- 0.66

- 0.31
- 50
- 0

- ω1
- ω2
- ω3

- P**k****k**

**k**

- P**k**slack**k**(s)

**k**

**s**

- P**k**H**k**

**k**

- - 6.46

- =

- ^

- + t n

## Case Study

**(6.77, 50, ****0****)**

**(****6.46****, 54, 0.09)**

**(11.5, ****0****, 0.662)**

- Pareto Surface

**Expected makespan**

**Satisfying demand**

**Robustness**

**Case Study**

**Notes:**
We use the same case study and the demand varies from 20 to 100. These 3 points are individual minima points and this is the pareto surface that covers 3 schedules. As we can see, the Pareto surface shows the tradeoff between reducing makespan, satisfying demand and better robustness. The decision maker can choose his solution according to this provided information.

## Case Study: Crude Oil Unloading

**Crude Oil **
**Marine Vessels**

**Storage**** **
**Tanks**

**Charging**
**Tanks**

**Crude**
**Distillation**
**Units**

**Other**
**Production**
**Units**

**Component**
**Stock Tanks**

**Blend**
**Header**

**Finished **
**Product**
**Tanks**

**Lifting/Shipping**
**Points**

- Problem 1: Crude-oil Unloading
- 	  and Mixing

- Problem 2: Production Stage

- Problem 3: Product Blending
- 	  and Distribution

## Case Study: Crude Oil Unloading

- (87.23, 0.06)

- (66.19, 2.14)

## Acknowledgements

- Financial Support
- BOC Gases, NSF CAREER Award (CTS-9983406), Petroleum Research Fund, Office of Naval Research

## Our Web Page:

**http://sol.rutgers.edu/staff/marianth**

**THANKS******