Distributionally Robust Optimization for the Newsvendor Problem
This example demonstrates how to solve the newsvendor problem based on the distributionally robust optimization framework (Bertsimas, Sim and Zhang 2014) with the XProg toolbox. Let \(\tilde{\pmb{z}}\in\mathbb{R}^N\) denote the vector of random variables, and the probability distribution is characterized the following ambiguity set \(\mathbb{F}\).
\begin{align}
\mathbb{F}=\left\{
\mathbb{P}\in\mathcal{P}_0\left(\mathbb{R}^N\right):
\begin{array}
\tilde{\pmb{z}}\in\mathbb{R}^N \\
\mathbb{E}_{\mathbb{P}}\left\{\tilde{\pmb{z}}\right\}=\pmb{\mu} \\
\mathbb{E}_{\mathbb{P}}\left\{\tilde{z}_i^2\right\}\leq \mu_i^2+\sigma_i^2,~~~~~\forall i=1,2,...,N \\
\mathbb{E}_{\mathbb{P}}\left\{\sum\limits_{i=1}^Np_i\cdot\max\left\{\mu_i-\tilde{z}_i,~0\right\}\right\}\leq \psi \\
\mathbb{P}\left\{\pmb{0}\leq \tilde{\pmb{z}}\leq \bar{\pmb{z}}\right\}=1
\end{array}
\right\}
\end{align} The extended ambiguity set is then formulated as below.
\begin{align}
\mathbb{G}=\left\{
\mathbb{Q}\in\mathcal{P}_0(\mathbb{R}^N\times\mathbb{R}^{N+1}\times\mathbb{R}^N):
\begin{array}
\tilde{\pmb{z}}\in\mathbb{R}^N,~\tilde{\pmb{u}}\in\mathbb{R}^{N+1},~\tilde{\pmb{v}}\in\mathbb{R}^N\\
\mathbb{E}_{\mathbb{Q}}\left\{\tilde{\pmb{z}}\right\}=\pmb{\mu} \\
\mathbb{E}_{\mathbb{Q}}\left\{\tilde{u}_i\leq \mu_i^2+\sigma_i^2\right\},~~~~~\forall i=1,2,...,N \\
\mathbb{E}_{\mathbb{Q}}\left\{\tilde{u}_{N+1}\right\}\leq \psi \\
\mathbb{P}\left\{\left(\tilde{\pmb{z}}, \tilde{\pmb{u}}, \tilde{\pmb{v}}\right)\in\bar{\mathcal{Z}}\right\}=1 \\
\end{array}
\right\}
\end{align} where the extended support set \(\mathcal{Z}\) is expressed as the following equation.
\begin{align}
\bar{\mathcal{Z}}=\left\{
\left(\pmb{z},\pmb{u},\pmb{v}\right)\in\mathbb{R}^N\times\mathbb{R}^{N+1}\times\mathbb{R}^N:
\begin{array}
\pmb{0}\leq \pmb{z}\leq \bar{\pmb{z}}\\
z_i^2\leq u_i, ~~~~~~\forall i=1,2,...,N\\
\pmb{p}^T\pmb{u}\leq u_{N+1} \\
\pmb{\mu}-\pmb{z}\leq \pmb{v} \\
\pmb{v}\geq \pmb{0} \\
\end{array}
\right\}
\end{align} Let \(\bar{\pmb{y}}(\pmb{z},\pmb{u},\pmb{v})\) be the decision rule that approximates the second-stage recourse decisions, which can expressed as the linear affine function below.
\begin{align}
\bar{\pmb{y}}(\pmb{z},\pmb{u},\pmb{v})=\pmb{y}^0+\sum\limits_{i=1}^N\pmb{y}_i^zz_i+\sum\limits_{j=1}^{N+1}\pmb{y}_j^uu_j+\sum\limits_{k=1}^{K}\pmb{y}_i^vv_k
\end{align} By applying the linear decision rule approximation, the newsvendor problem can be formulated as follows.
\begin{align}
\min~&-\pmb{p}^T\pmb{x}+\max\limits_{\mathbb{Q}\in\mathbb{G}}\mathbb{E}_{\mathbb{Q}}\left\{\pmb{p}^T\bar{\pmb{y}}\left(\tilde{\pmb{z}}, \tilde{\pmb{u}}, \tilde{\pmb{v}}\right)\right\} \\
\text{s.t.}~&\pmb{c}^T\pmb{x}\leq \Gamma\\
&\pmb{x}\geq \pmb{0} \\
&\bar{\mathbb{y}}(\pmb{z},\pmb{u},\pmb{v})\geq \pmb{x}-\pmb{z}, ~~~~~\forall \left(\pmb{z}, \pmb{u}, \pmb{v}\right)\in\bar{\mathcal{Z}} \\
&\bar{\mathbb{y}}(\pmb{z},\pmb{u},\pmb{v})\geq \pmb{0}, ~~~~~~~~~~~~\forall \left(\pmb{z}, \pmb{u}, \pmb{v}\right)\in\bar{\mathcal{Z}}
\end{align} where \(\pmb{x}\in\mathbb{R}^N\) is the vector of the first-stage decisions, The budget constant \(\Gamma\) is set to be \(500\), and another constant (\psi=100\). The data of the price coefficient vector \(\pmb{p}\), the cost coefficient vector \(\pmb{c}\), and the upper bounds of random varaibles, can be found in reference (Bertsimas, Sim and Zhang 2014). The corresponding code is given below, and the objective value
This example demonstrates how to solve the newsvendor problem based on the distributionally robust optimization framework (Bertsimas, Sim and Zhang 2014) with the XProg toolbox. Let \(\tilde{\pmb{z}}\in\mathbb{R}^N\) denote the vector of random variables, and the probability distribution is characterized the following ambiguity set \(\mathbb{F}\).
\begin{align}
\mathbb{F}=\left\{
\mathbb{P}\in\mathcal{P}_0\left(\mathbb{R}^N\right):
\begin{array}
\tilde{\pmb{z}}\in\mathbb{R}^N \\
\mathbb{E}_{\mathbb{P}}\left\{\tilde{\pmb{z}}\right\}=\pmb{\mu} \\
\mathbb{E}_{\mathbb{P}}\left\{\tilde{z}_i^2\right\}\leq \mu_i^2+\sigma_i^2,~~~~~\forall i=1,2,...,N \\
\mathbb{E}_{\mathbb{P}}\left\{\sum\limits_{i=1}^Np_i\cdot\max\left\{\mu_i-\tilde{z}_i,~0\right\}\right\}\leq \psi \\
\mathbb{P}\left\{\pmb{0}\leq \tilde{\pmb{z}}\leq \bar{\pmb{z}}\right\}=1
\end{array}
\right\}
\end{align} The extended ambiguity set is then formulated as below.
\begin{align}
\mathbb{G}=\left\{
\mathbb{Q}\in\mathcal{P}_0(\mathbb{R}^N\times\mathbb{R}^{N+1}\times\mathbb{R}^N):
\begin{array}
\tilde{\pmb{z}}\in\mathbb{R}^N,~\tilde{\pmb{u}}\in\mathbb{R}^{N+1},~\tilde{\pmb{v}}\in\mathbb{R}^N\\
\mathbb{E}_{\mathbb{Q}}\left\{\tilde{\pmb{z}}\right\}=\pmb{\mu} \\
\mathbb{E}_{\mathbb{Q}}\left\{\tilde{u}_i\leq \mu_i^2+\sigma_i^2\right\},~~~~~\forall i=1,2,...,N \\
\mathbb{E}_{\mathbb{Q}}\left\{\tilde{u}_{N+1}\right\}\leq \psi \\
\mathbb{P}\left\{\left(\tilde{\pmb{z}}, \tilde{\pmb{u}}, \tilde{\pmb{v}}\right)\in\bar{\mathcal{Z}}\right\}=1 \\
\end{array}
\right\}
\end{align} where the extended support set \(\mathcal{Z}\) is expressed as the following equation.
\begin{align}
\bar{\mathcal{Z}}=\left\{
\left(\pmb{z},\pmb{u},\pmb{v}\right)\in\mathbb{R}^N\times\mathbb{R}^{N+1}\times\mathbb{R}^N:
\begin{array}
\pmb{0}\leq \pmb{z}\leq \bar{\pmb{z}}\\
z_i^2\leq u_i, ~~~~~~\forall i=1,2,...,N\\
\pmb{p}^T\pmb{u}\leq u_{N+1} \\
\pmb{\mu}-\pmb{z}\leq \pmb{v} \\
\pmb{v}\geq \pmb{0} \\
\end{array}
\right\}
\end{align} Let \(\bar{\pmb{y}}(\pmb{z},\pmb{u},\pmb{v})\) be the decision rule that approximates the second-stage recourse decisions, which can expressed as the linear affine function below.
\begin{align}
\bar{\pmb{y}}(\pmb{z},\pmb{u},\pmb{v})=\pmb{y}^0+\sum\limits_{i=1}^N\pmb{y}_i^zz_i+\sum\limits_{j=1}^{N+1}\pmb{y}_j^uu_j+\sum\limits_{k=1}^{K}\pmb{y}_i^vv_k
\end{align} By applying the linear decision rule approximation, the newsvendor problem can be formulated as follows.
\begin{align}
\min~&-\pmb{p}^T\pmb{x}+\max\limits_{\mathbb{Q}\in\mathbb{G}}\mathbb{E}_{\mathbb{Q}}\left\{\pmb{p}^T\bar{\pmb{y}}\left(\tilde{\pmb{z}}, \tilde{\pmb{u}}, \tilde{\pmb{v}}\right)\right\} \\
\text{s.t.}~&\pmb{c}^T\pmb{x}\leq \Gamma\\
&\pmb{x}\geq \pmb{0} \\
&\bar{\mathbb{y}}(\pmb{z},\pmb{u},\pmb{v})\geq \pmb{x}-\pmb{z}, ~~~~~\forall \left(\pmb{z}, \pmb{u}, \pmb{v}\right)\in\bar{\mathcal{Z}} \\
&\bar{\mathbb{y}}(\pmb{z},\pmb{u},\pmb{v})\geq \pmb{0}, ~~~~~~~~~~~~\forall \left(\pmb{z}, \pmb{u}, \pmb{v}\right)\in\bar{\mathcal{Z}}
\end{align} where \(\pmb{x}\in\mathbb{R}^N\) is the vector of the first-stage decisions, The budget constant \(\Gamma\) is set to be \(500\), and another constant (\psi=100\). The data of the price coefficient vector \(\pmb{p}\), the cost coefficient vector \(\pmb{c}\), and the upper bounds of random varaibles, can be found in reference (Bertsimas, Sim and Zhang 2014). The corresponding code is given below, and the objective value
N=10; Price=[10.00 11.00 11.41 11.73 12.00 12.24 12.45 12.65 12.83 13.00]'; % vector of price coefficients: p Cost =[2.00 2.71 3.00 3.23 3.41 3.58 3.73 3.87 4.00 4.12 ]'; % vector of cost coefficients: c mu =[30.00 35.00 40.00 45.00 50.00 55.00 60.00 65.00 70.00 75.00]'; % mean values of random variables z sig =(30:-1.5:16.5)'; % standard deviations of random variables z barZ =100*ones(10,1); % upper bounds of random variables z Psi =100; Gamma=500; model=xprog('newsvendor'); % create a model, named 'newsvendor' x=model.decision(N); % here-and-now decisions x y=model.recourse(N); % decision rule y for recourse decisions z=model.random(N); % define random variables z u=model.random(N+1); % define auxiliary variables u v=model.random(N); % define auxiliary variables z y.depend(z); % define dependency of y on z y.depend(u); % define dependency of y on u y.depend(v); % define dependency of y on v % extended ambiguity set is defined below model.uncertain(expect(z)==mu); % the 2nd line of set G mu2 =mu.*mu; % mu^2 sig2=sig.*sig; % sigma^2 model.uncertain(expect(u(1:N)')<=mu2+sig2); % the 3rd line of set G model.uncertain(expect(u(N+1))<=Psi); % the 4th line of set G model.uncertain(z>=0); % lower bound of z model.uncertain(z<=barZ); % upper bound of z model.uncertain(z.^2<=u(1:N)'); % the 2nd line of set \bar{Z} model.uncertain(u(N+1)>=Price'*v); % the 3rd line of set \bar{Z} model.uncertain(v>=mu-z); % the 4th line of set \bar{Z} model.uncertain(v>=0); % the 5th line of set \bar{Z} model.max(Price'*x-expect(Price'*y)); % define objective function model.add(Cost'*x<=Gamma); % the 1st constraint: budget model.add(x>=0); % the 2nd constraint: nonnegtivity model.add(y>=x-z); % the 3rd constraint model.add(y>=0); % the 4th constraint model.solve; % solve the problem
The solution of this problem is provided in the table below.
Objective | x1 | x2 | x3 | x4 | x5 | x6 | x7 | x8 | x9 | x10 |
1850.10 | 30.00 | 35.00 | 40.00 | 45.00 | 23.40 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 |