610 likes | 897 Views
Short-Term Scheduling under Uncertainty. Marianthi Ierapetritou Department Chemical and Biochemical Engineering Piscataway, NJ 08854-8058. Process Operations Decision Making. Uncertainty Complexity. Short-term Scheduling. Time Horizon. Production Planning. Supply Chain Management.
E N D
Short-Term Scheduling under Uncertainty Marianthi Ierapetritou Department Chemical and Biochemical Engineering Piscataway, NJ 08854-8058
Process Operations Decision Making Uncertainty Complexity Short-term Scheduling Time Horizon Production Planning Supply Chain Management Online Control Objective • Identify and reduce bottlenecks at different levels • Integration of the whole decision-making process Opportunity for Optimization
Uncertain Parameters 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
Short-term Scheduling Process Plant Optimal Schedule Given: Determine: Raw Materials, Required Products, Task Sequence, Production Recipe, Unit Capacity Exact Amounts of material Processed Scheduling objectives : EconomicMaximize Profit, Minimize Operating Costs, Minimize Inventory Costs Time BasedMinimize Makespan, Minimize Tardiness
Continuous Time Formulation Discrete time formulation T1 T1 T2 T2 T1 T2 T2 T2 “Real” time schedule Continuous time formulation T1 T1 T2 T2 Event Points: When a tasks begins • Binary variables to allocate tasks to resources • Continuous variables to represent timing and material variables • Mixed IntegerLinear Programming Models • Smaller models that are computationally efficient and tractable
Deterministic Scheduling Formulation minimize H or maximizeprice(s)d(s,n) subject towv(i,j,n) 1 st(s,n) = st(s,n-1) – d(s,n) + Pb(i,j,n-1) + cb(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 Objective Function s Allocation Constraints (i,j) Material Balances Capacity Constraints Demand Constraints n Duration Constraints M.G.Ierapetritou and C.A.Floudas. Effective continuous-time formulation for short-term scheduling. 1. Multipurpose batch processes. 1998
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 Deterministic Schedule 55.56 separation 44.44 74.07 50.93 reaction 1 reaction 1 reaction 3 demand (product 2) = 50 74.07 4.63 50.93 reaction 2 reaction 3 reaction 2 E(makespan) = 8.15hr 50.00 heating Standard Deviation = 2.63 0 1 2 3 4 5 6 7 8 Robust Schedule 55.56 separation 64.60 10.03 50.93 reaction 3 reaction 1 reaction 2 demand (product 2) = 50*(1 + 60%) 10.40 64.04 4.63 50.93 reaction 1 reaction 2 reaction 3 reaction 2 E(makespan) = 7.24hr 50.00 heating Standard Deviation = 0.29 0 1 2 3 4 5 6 7 8
Literature Review: Representative Publications • Reactive Scheduling Handles uncertainty by adjusting a schedule upon realization of the uncertain parameters or occurrence of unexpected events • 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 • 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 • Fuzzy Programming Considers random parameters as fuzzy numbers and the constraints are treated as fuzzy sets • 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 • Robust Optimization Produces “robust” solutions that are immune against uncertainties • 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 • MILP Sensitivity Analysis Utilizes MILP sensitivity analysis methods to investigate the effects of uncertain parameters and provide a set of alternative schedules • Z.Jia and M.G.Ierapetritou. Short-term Scheduling under Uncertainty Using MILP Sensitivity Analysis. Ind. Eng. Chem. Res. 2004, 43, 3782
Disruptive Events Rush Order Arrivals Order Cancellations Machine Breakdowns Not much information is available REACTIVE SCHEDULING Parameter Uncertainty Processing times Demand of products Prices Information is available PREVENTIVE SCHEDULING Uncertainty in Scheduling
Preventive Scheduling New alternative schedules MILP sensitivity analysis framework Data perturbation LB/UB on objective function Deterministic schedule Robust optimization method A set of solutions represent trade-off between various objectives model robustness solution robustness
Preventive Scheduling • MILP Sensitivity Analysis minimize H or maximizeprice(s)d(s,n) subject towv(i,j,n) 1 st(s,n) = st(s,n-1) – d(s,n) + Pb(i,j,n-1) + cb(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 • Robust Optimization
Questions to Address • What is the effect of processing time at the objective value? 55 15 mixing mixing 55 15 reaction reaction 50 20 separation separation 0 2 4 8 10 6 H (time horizon) • Can the schedule accommodate the demand fluctuation? • How the capacity of the units affect the production objective?
j = 1,…,n xj {ujP,…,ujP} Inference-based MILP Sensitivity Analysis 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 Bound z z*- z holdsif there ares1P,…,snPthat satisfy: - for the perturbations A and a - for the perturbations c iP AijujP + sj(uj – uj) - iai rP sjPiPAij, sjP -qjP, j = 1,…,n rP = -qjPujP +Pa – zP +zP cjujP - sjP(ujP – ujP) -rP sjP -cj,sjP-qjP, j = 1,…,n qjP = iPAij - iPcj *M.W.Dawande and J.N.Hooker, 2000
Solve relaxed LP at the leaf nodes with perturbed data Identify the feasible schedules by examining the B&B tree • Robustness • Nominal performance • Average performance Evaluate the alternative schedules Proposed Uncertainty Analysis Approach Solve the deterministic scheduling problem using B&B tree MILP Sensitivity Analysis Extract information from the leaf nodes • Range of objectivechange for certain parameter change
Robustness Estimation Makespan minimization is considered as the objective Obtain sequence of tasks from original schedule Generate random demands in expected range Makespan to meet a particular demand is found using the sequence of tasks derived from original schedule Binary variablescorresponding to allocation of tasks arefixed Batch sizesandStartingandFinishingtimes of tasks are allowed to vary
Inventory of raw materials intermediates etc. ROLLOVER Hmax Hinf Meets unsatisfied demand Meets maximum possible demand Total makespan Hcorr = Hmax + Hinf Robustness under Infeasibility Corrected Standard Deviation: Hact = Hp if scenario is feasible = Hcorr if the scenario is infeasible J.P.Vin and M.G.Ierapetritou. Robust short-term scheduling of multiproduct batch plants under demand uncertainty. 2001
Case Study 1 S1 S2 S3 S4 mixing reaction purification Effect of demand d~[20, 100] -0.097 d H dnom = 50 Hnom = 9.83h d’= 80 H’ Hnom + 0.097d =12.73h 3.0 B&B tree with nominal demand wv(i1,j1,n0) 1 0 5.17 3.0 wv(i1,j1,n1) 1 0 0 1 7.65 5.83 7.16 infeasible wv(i2,j2,n1) 1 0 1 0 1 0 8.14 7.16 7.16 10.16 8.33 infeasible wv(i2,j2,n2) 1 0 1 0 0 1 0 1 9.87 8.83 9.98 8.83 9.83 9.83 infeasible infeasible wv(i3,j3,n2) (Schedule 1)
Case Study 1 1 0 1 0 1 0 1 0 (Schedule 2) (Schedule 3) 3.0 wv(i1,j1,n0) 1 0 5.17 3.0 wv(i1,j1,n1) 1 0 0 1 7.65 5.83 7.16 infeasible wv(i2,j2,n1) 1 0 1 0 1 0 8.14 7.16 7.16 10.16 8.33 infeasible wv(i2,j2,n2) 1 0 1 0 0 1 0 1 9.87 8.83 9.98 8.83 9.83 9.83 (12.13) (17.97) (12.73) infeasible infeasible wv(i3,j3,n2) (Schedule 1) wv(i3,j3,n3) schedule 2 schedule 3 schedule 1 Schedule Evaluation 10.77 10.91 Hnom(h) 9.83 11.56 Havg(h) 14.20 11.79 SDcorr 1.61 5.52 2.17
Case Study 1 schedule 3 schedule 2 (optimal when d ≥ 50) schedule 1 (optimal when d ≤ 50)
1 1 1 1 0 1 0 1 0 (Schedule 2) (Schedule 3) Case Study 1 Effect of processing timeT(i1,j1) ~ [2.0, 4.0] Tnom = 3.0 profitnom = 71.52 profit’ profitnom + 24.48T = 47.04 T’= 4.0 schedule 2 schedule 3 schedule 1 100 wv(i1,j1,n0) 65.27 65.27 profitnom 71.52 1 0 100 50 profitavg 66.98 65.17 64.61 wv(i1,j1,n1) 1 0 SDcorr 10.49 9.33 26.9 100 100 wv(i2,j2,n1) 1 0 1 0 100 50 96.05 50 wv(i2,j2,n2) 1 0 1 0 75 78.42 72.46 75 (75) (62.11) (75) wv(i3,j3,n2) 1 0 78.42 50 wv(i3,j3,n3) 1 0 71.52 50 (Schedule 1)
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 z() = min cTx + dTy subject to Ax + Dy b xL x xU L U x Rm, y (0,1)T b = b0 + r b[b0+Lr, b0+Ur] solved at b = b0+Lr optimal solution (x*,y*) Fix integer variables at y* LP sensitivity analysis: 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 z() = z0 + L L’ U’ U Integer cut to exclude current optimal solution iF1 iF 0 break point ’, new optimal solution (x*’, y*’) A.Pertsinidis et al. Parametric optimization of MILP programs and a framework for the parametric optimization of MINLPs. 1998
Multiparametric MILP Solve the fully relaxed problem Select a branching variable Based on simplex algorithm, check the neighboring bases of the LP tableau Solve the mpLP at the nodes Compare the solution with the current UB, update the optimal function in the uncertain space J. Acevedo and E.N.Pistikopoulos. A Multiparametric Programming Approach for Linear Process Engineering Problems under Uncertainty. 1997
Shortcomings of Existing Approach • Solve mpLP at every node in the B&B tree during the branch and bound procedure: can be a computationally expensive effort • mpLP approach requires retrieving the LP tableaus and visiting the neighbor bases
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 Solve the original problem at the nominal value using a branch and bound method Find Δamax that leaves the structure of the tree unchanged Collect zp, λp at each leaf node p For a = Δamax + ε update the B&B tree
zp – z0 case 3 λp – λ0 case 1 * * case 2 Determine Δamax and Update the B&B Tree where node 0 is the optimal node Δamax = min{Δabasis, min{ }} P Update the B&B tree at a = Δamax +ε 0 The new optimal node can be: 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
Proposed Approach: Multiple Uncertain Parameters Solve the original problem at the nominal value using a branch and bound method mpLP at the leaf nodes mpLP algorithm Compare the critical regions with the current upper bounds Update the B&B tree
mpLP Algorithm at the Leaf Nodes At each iteration, solve maximize cx*– z subject to z = max{z(k) + λ(k)θa + β(k)θb} a0 ≤ θa≤ a0 + Δa b0≤ θb≤ b0 + Δb Bilevel Linear Programming where x* = argmin cx subject to Ax ≥θ cx* cx* current optimal functions max{z(k) + λ(k)θ} θ
mpLP Algorithm at the Leaf Nodes check if there is any point at which max{z(k) + λ(k)θa + β(k)θb} is less than cx* 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 Include an additional constraint z z(k+1) + λ(k+1)θa + β(k+1)θb to above problem Stop when the objective = 0
Solve the Bilevel Programming Problem Convert the relaxed LP to its dual form mincx s.t. Ax θ a0 ≤ θa≤ a0 + Δa b0≤ θb≤ b0 + Δb x 0 maxθy s.t. ATy ≤c a0 ≤ θa≤ a0 + Δa b0≤ θb≤ b0 + Δb Y 0 Strong Duality Theorem Replace the inside optimization problem by its dual Solve with global optimization solver BARON Bilinear objective Linear constraints 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
Compare the Optimal Functions of the Leaf Nodes 1 2 CR2(1) CR2(2) CR1(1) CR1UB CR2(3) z1UB = z1*UB + λ1UBθa + β1UBθb z2(2)= z2*(2) + λ2(2)θa + β2(2)θb CR1UBCR2(2)= CRint
Compare the Optimal Functions of the Leaf Nodes Compare optimal function z1UB and z2(2) in region CRint 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. 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 • 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 • The new mpLP approach can efficiently determine the optimal functions and critical regions without retrieving the LP tableaus and visiting the neighbor bases
Case Study: Single Uncertain Parameter minz = 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: Δa = 0 x3 = 1 x3 = 0 2 12 x4 = 1 x4 = 0 0 3 11.5* 12 (0,1,1) Step 2: Linear sensitivity analysis on node 0 -- Δamax = 0 For Δa = Δamax + ε, node 0 yields noninteger solution. Update the B&B tree. Step 3:
zp – z0 λp – λ0 0.3 0.2 0.2 0.2 , , 2 3 3 3 Case Study: Single Uncertain Parameter Step 2: x3 = 1 x3 = 0 2 Δamax = min{Δabasis, min{ }} 12 x4 = 1 x4 = 0 P λ = 0 = min{ 2, min { }} 0 3 12 x5 = 1 x5 = 0 λ = 0 = 11.8* 12.1 λ = 3 λ = 1 Step 3: For Δa = Δamax + ε, node 0 is intersected by node 2 & 3. Update the B&B tree x3 = 1 x3 = 0 2 1 0 1 0 ... Step 2: 12.2 4 12* 0 3 λ = 1 λ = 2 ... 1 0 1 0 Step 3: 5 12.1 12.2 12*
Case Study: Single Uncertain Parameter Final optimal solution: 14 13.5 (1,0,1) 13 z* 12.5 (1,0,1) (1,0,1) & (0,0,0) 12 (0,1,1) 11.5 0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2 Δa
Case Study: Multiple Uncertain Parameters minz = -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, θ2) = (0, 0) y1 = 0 y1 = 1 y2 = 0 y2 = 1 2 1 -8 -70.5*
Case Study: Multiple Uncertain Parameters mpLP on the leaf nodes Step 2: z1(1) = -70.5 – 4.3333θ1 – 0.1667θ2 at (θ1, θ2) = (0, 0) Node 1 max(-13-θ1)d1 – 20d2 + (-121-θ2)d3 + 8d4 – d7 – d8 - z 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 obj = 16.74 (θ1, θ2) = (10, 0) z1(2) = -97.0909 – 0.1667θ2 at (θ1, θ2) = (10, 0)
Case Study: Multiple Uncertain Parameters max(-13-θ1)d1 – 20d2 + (-121-θ2)d3 + 8d4 – d7 – d8 - z 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 The two critical regions are divided by z1(1) ≥ z1(2) 0.07333θ1 – 0.00333θ2 ≤ 0.45
Case Study: Multiple Uncertain Parameters z1(2) = -97.0909 – 0.3636θ2 z1(1) = -70.5 – 4.3333θ1 – 0.1667θ2 0.07333θ1 – 0.00333θ2 ≤ 0.45 0.07333θ1 – 0.00333θ2 ≥ 0.45 CR1(1)= CR1(2)= θ2 ≤ 10 θ1, θ2 ≤ 10 z2 = 8 Node 2 θ1 ≤ 10 CR2 = θ2 ≤ 10 Step 3: Compare critical regions and determine optimal functions CR1(1)n CR2 = CR1(1) CR1(1) and CR2
Case Study: Multiple Uncertain Parameters minЄ s.t. -70.5 – 4.3333θ1 - 0.1667θ2 + Є =-8 0.07333θ1 – 0.00333θ2 ≤ 0.45 θ2 ≤ 10 CR1(1) z1(1)≤ z2 Є < 0, CR1(2)n CR2 = CR1(2) CR1(2) and CR2 minЄ s.t. -97.0909 – 0.3636θ2 + Є =-8 0.07333θ1 – 0.00333θ2 ≥ 0.45 0 ≤ θ1, θ2 ≤ 10 CR1(2) z1(2)≤ z2 is redundant in CR1(2) Є > 0,
Case Study: Multiple Uncertain Parameters Final optimal solution: y* = (1, 1) z1(θ) = -70.5 – 4.3333θ1 – 0.1667θ2 0.07333θ1 – 0.00333θ2 ≤ 0.45 CR1(1)= θ2 ≤ 10 z2(θ) = -97.0909 – 0.3636θ2 0.07333θ1 – 0.00333θ2 ≥ 0.45 CR1(2)= θ1, θ2 ≤ 10
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 Unlike the case of uncertain RHS, where the optimal objective value z* = max{z(k) + λ(k)θa + β(k)θb} Here, z* = min{z(k) + λ(k)c1 + β(k)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 ≥θ maximize z –(c + Δc)x subject to Ax ≥ θ z ≤ z(k) + λ(k)c1 + β(k)c2 c10 ≤ c1≤ c10 + Δc1 c20≤ c2≤ c20 + Δc2 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 At mpLP procedure, need to solve 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 ≥θ
Solve Bilevel Linear Programming Problem The most popular method is “Kuhn-Tucker” approach min F(x, y) = c1x + d1y subject to A1x + B1y ≤ b1 xЄX min f(x, y) = c2x + d2y subject to A2x + B2y ≤ b2 Replace with its KKT condition and add to the upper level problem yЄY