Robust Price Optimization of Multiple Products under Interval Uncertainties

For our MPPO, we assume a linear price-response function as suggested by Phillips [23]:

$$d(p)=\alpha+\beta p$$

(2.1)

where $\alpha>0$ and $\beta<0$ are the respective intercept and slope of the linear price-response function $d(p)$ which represents the sales volume as a function of price. The satiating price $P$, defined as the maximum allowed price at which demand drops to zero, is given by $P=-\alpha/\beta$. Phillips [23] points out that the elasticity of the linear price-response function is -$\beta p/(\alpha+\beta p)$ which ranges from 0 at $p=0$ and approaches infinity as $p$ approaches $P$.

In reality, it is well known that the demand of a product is often influenced by not only the price of itself, but also the prices of its complementary and substitute products. Therefore we extend the linear function (2.1) to incorporate the prices of complementary and substitute products. Let $I$ represent the set of product indexes. The extended demand function for product $i$ can be expressed by introducing additional price variables along with complementary/substitute products:

$$d_{i}(p)=\alpha_{i}+\beta_{i}p_{i}+\sum_{j\in C_{i}}{\gamma_{ij}p_{j}}$$

(2.2)

where $p_{j}$ is the price of the complementary/substitute product $j$, $\gamma_{ij}$ represents the corresponding cross-effect of product $j$ on own product $i$, and $C_{i}$ is the set of associated complementary/substitute products for product $i$.

The MPPO is to maximize the total revenue subject to a set of constraints:

$$\begin{array}[]{rl}\max\limits_{p}&\sum\limits_{i\in I}{p_{i}\cdot(\alpha_{i}+\beta_{i}p_{i}+\sum\limits_{j\in C_{i}}{\gamma_{ij}p_{j}})\vspace{1mm}}\\ \mathop{\rm s.t.}&Ap\leq b,\vspace{1mm}\\ &l\leq p\leq u,\end{array}$$

(2.3)

where $p=(p_{1},p_{2},\ldots,p_{n})\in\mathbb{R}^{n}$ (hence $n=|I|$) , $A\in\mathbb{R}^{m\times n}$, $b\in\mathbb{R}^{m}$, and $l,u\in\mathbb{R}_{+}^{n}$ are finite. The constraints in (2.3) defined by matrix $A$ and vector $b$ represent business rules in practice.

We assume that the above problem is well-defined and feasible. We also assume that for $i\in I$, $\beta_{i}<0$, so the demand of a product decreases when its price increases, and $\gamma_{ij}<0$ for any complementary product $j\in C_{i}$ and $\gamma_{ij}>0$ for any substitute product $j\in C_{i}$.

Thus, the objective function

$$f(p;\beta)=\sum_{i\in I}{\left(\alpha_{i}p_{i}+\beta_{i}p_{i}^{2}+\sum_{j\in C_{i}}{\gamma_{ij}p_{i}p_{j}}\right)},$$

is, in general, a nonconvex quadratic function with bilinear terms.

For the quadratic and bilinear terms, $p_{i}^{2}$ and $p_{i}p_{j}$, in (2.3), we introduce two new sets of variables $x_{i}=-p_{i}^{2}$ and $y_{ij}=p_{i}\cdot p_{j}$, respectively. With this substitution, (2.3) is transformed into:

	$\displaystyle\max\limits_{p,x,y}\quad$	$\displaystyle\sum\limits_{i\in I}{\left(\alpha_{i}p_{i}-\beta_{i}x_{i}+\sum_{j\in C_{i}}{\gamma_{ij}y_{ij}}\right)}\vspace{1mm}$			(2.4a)
	$\displaystyle\mathop{\rm s.t.}\quad$	$\displaystyle Ap\leq b,$			(2.4b)
		$\displaystyle x_{i}\leq-p_{i}^{2},$	$\displaystyle i\in I,$		(2.4c)
		$\displaystyle y_{ij}=p_{i}p_{j},$	$\displaystyle i\in I,j\in C_{i},$		(2.4d)
		$\displaystyle l\leq p\leq u.$			(2.4e)

where the objective function (2.4a) is now linear, but the constraints (2.4c) and (2.4d) contain quadratic and bilinear terms, respectively.

Notice that for a given $i\in I$, since $\beta_{i}<0$ and $p_{i}>0$, $x_{i}<0$ and $-\beta_{i}x_{i}<0$. Because we are dealing with a maximization problem, it follows that $x_{i}$ will attain its least negative value at $x_{i}=-p^{2}_{i}$ at optimality. Therefore, Problem (2.4) and Problem (2.3) are equivalent and consequently, we can deal with Problem (2.4) from this point on.

We now apply relaxations to Problem (2.4) to address the nonlinear constraints (2.4c) and (2.4d) and transform the problem into a more tractable linear program. First, we replace the quadratic term $-p_{i}^{2}$ by its piecewise linear approximation by discretizing its domain $p_{i}$ into $r$ segments:

$$l_{i}=\bar{p}_{i0}\leq\bar{p}_{i1}\leq\ldots\leq\bar{p}_{ir}=u_{i}.$$

Then at the $k$th knot:

$$\begin{array}[]{rcl}p_{i}^{2}&\approx&-\bar{p}_{ik}^{2}+(-2\bar{p}_{ik})(p_{i}-\bar{p}_{ik})\\ &=&\bar{p}_{ik}^{2}-2\bar{p}_{ik}p_{i}.\end{array}$$

Therefore, we can relax constraint (2.4c) with the following set of $r$ supporting hyperplanes of $p_{i}^{2}$:

$$x_{i}\leq\bar{p}_{ik}^{2}-2\bar{p}_{ik}p_{i}\quad i\in C,k\in\{1,\ldots,r\}.$$

Figure 1 illustrates how this linearization technique to approximate the quadratic function $x_{i}=-p_{i}^{2}$.

Next, for constraint (2.4d), we relax the bilinear terms by introducing McCormick envelopes ([19, 2]) as follows:

	$\displaystyle y_{ij}\geq l_{j}p_{i}+l_{i}p_{j}-l_{i}l_{j},$
	$\displaystyle y_{ij}\leq u_{j}p_{i}+l_{i}p_{j}-l_{i}u_{j},$
	$\displaystyle y_{ij}\leq l_{j}p_{i}+u_{i}p_{i}-l_{j}u_{i},$
	$\displaystyle y_{ij}\geq u_{j}p_{i}+u_{i}p_{j}-u_{i}u_{j}.$

With the above relaxations, Problem (2.4) becomes:

	$\displaystyle\max\limits_{p,x,y}\quad$	$\displaystyle\sum\limits_{i\in I}{\left(\alpha_{i}p_{i}-\beta_{i}x_{i}+\sum_{j\in C_{i}}{\gamma_{ij}y_{ij}}\right)}\vspace{1mm}$			(2.5a)
	$\displaystyle\mathop{\rm s.t.}\quad$	$\displaystyle Ap\leq b,$			(2.5b)
		$\displaystyle x_{i}\leq\bar{p}_{ik}^{2}-2\bar{p}_{ik}p_{i}$	$\displaystyle i\in I,k\in\{1,\ldots,r\},$		(2.5c)
		$\displaystyle y_{ij}\geq l_{j}p_{i}+l_{i}p_{j}-l_{i}l_{j},$	$\displaystyle i\in I,j\in C_{i}$		(2.5d)
		$\displaystyle y_{ij}\leq u_{j}p_{i}+l_{i}p_{j}-l_{i}u_{j},$	$\displaystyle i\in I,j\in C_{i}$		(2.5e)
		$\displaystyle y_{ij}\leq l_{j}p_{i}+u_{i}p_{i}-l_{j}u_{i},$	$\displaystyle i\in I,j\in C_{i}$		(2.5f)
		$\displaystyle y_{ij}\geq u_{j}p_{i}+u_{i}p_{j}-u_{i}u_{j},$	$\displaystyle i\in I,j\in C_{i}$		(2.5g)
		$\displaystyle l\leq p\leq u.$			(2.5h)

The RO model is constructed by considering the worst-case scenario under all possible realizations of the uncertain parameter. Therefore, the RO counterpart of the objective function of Problem (2.5) is written as a Max-Min problem as follows:

		$\displaystyle\max\limits_{p,x,y}\quad$	$\displaystyle\min_{(\beta,\gamma)\in\mathcal{U}}\sum\limits_{i\in I}{\left(\alpha_{i}p_{i}-\beta_{i}x_{i}+\sum_{j\in C_{i}}{\gamma_{ij}y_{ij}}\right)}$
	$\displaystyle=$	$\displaystyle\max\limits_{p,x,y}\quad$	$\displaystyle\left(\sum\limits_{i\in I}{\alpha_{i}p_{i}}+\min_{(\beta,\gamma)\in\mathcal{U}}\sum\limits_{i\in I}{\Bigl{(}-\beta_{i}x_{i}+\sum_{j\in C_{i}}{\gamma_{ij}y_{ij}}\Bigr{)}}\right)$

Hence, Problem (2.5) under uncertainty becomes:

$$\begin{array}[]{rl}\max\limits_{p,x,y}&\left(\sum\limits_{i\in I}{\alpha_{i}p_{i}}+\min\limits_{(\beta,\gamma)\in\mathcal{U}}\sum\limits_{i\in I}{\Bigl{(}-\beta_{i}x_{i}+\sum\limits_{j\in C_{i}}{\gamma_{ij}y_{ij}}\Bigr{)}}\right)\vspace{3mm}\\ \mathop{\rm s.t.}&\text{\eqref{prob:reform2:const1} - \eqref{prob:reform2:const4}},\end{array}$$

(3.1)

where $\mathcal{U}$ is the uncertainty set of $(\beta,\gamma)$. We will describe a method for constructing the uncertainty set later.

By introducing an ancillary variable $\eta$, we reformulate the model (3.1) as

	$\displaystyle\max\limits_{p,x,y}\quad$	$\displaystyle\eta+\sum\limits_{i\in I}{\alpha_{i}p_{i}}$
	$\displaystyle\mathop{\rm s.t.}\quad$	$\displaystyle\eta\leq\min\limits_{(\beta,\gamma)\in\mathcal{U}}\left(\sum\limits_{i\in I}{\Bigl{(}-\beta_{i}x_{i}+\sum_{j\in C_{i}}{\gamma_{ij}y_{ij}}\Bigr{)}}\right)$		(3.2a)
		$\displaystyle\text{\eqref{prob:reform2:const1} - \eqref{prob:reform2:const4}}.$

We denote $S_{i}^{RO}$ and $T_{ij}^{RO}$ as two pre-determined numbers of uniformly distributed realizations $\hat{\beta}_{is}$ and $\hat{\gamma}_{ijs}$ for uncertain parameters $\beta_{i}$ and $\gamma_{ij}$, respectively, from their corresponding interval estimates. Given the set of realizations, we define the following attributes for the uncertainty set for each product $i$:

• •

  Mean of scenarios:

  $$\bar{\beta}_{i}=\frac{\sum\limits_{s=1}^{S_{i}^{RO}}{\hat{\beta}_{is}}}{S_{i}^{RO}}.$$

  $\bar{\beta}_{i}$ can be interpreted as the nominal value for the sensitivity parameter $\beta_{i}$. Likewise, for $\gamma_{ij}$, we have:

  $$\bar{\gamma}_{ij}=\frac{\sum\limits_{s=1}^{T_{ij}^{RO}}{\hat{\gamma}_{is}}}{T_{ij}^{RO}}.$$

• •

  Standard deviation of scenarios: The standard deviations of the samples of realizations for $\beta$ and $\gamma$ are calculated as follows.

  $$\tilde{\beta}_{i}=\sqrt{\frac{\sum\limits_{s=1}^{S_{i}^{RO}}{(\hat{\beta}_{is}-\bar{\beta}_{i})^{2}}}{S_{i}^{RO}-1}},$$

  and

  $$\tilde{\gamma}_{ij}=\sqrt{\frac{\sum\limits_{s=1}^{T_{ij}^{RO}}{(\hat{\gamma}_{ijs}-\bar{\gamma}_{ij})^{2}}}{T_{ij}^{RO}-1}}.$$

Now, for $k=\sum_{i\in I}{|C_{i}|}$ and writing $\boldsymbol{\beta}$ for the vectors of all $\beta$ and $\boldsymbol{\gamma}\in\mathbb{R}^{k}$ for the vector of all $\gamma_{ij}$, we construct the uncertainty set $\mathcal{U}_{\Delta}$ as follows:

	$\displaystyle\mathcal{U}_{\Delta}=\Bigl{\{}\left[\begin{array}[]{c}\boldsymbol{\beta}\\ \boldsymbol{\gamma}\end{array}\right]\in$	$\displaystyle\mathbb{R}^{n+k}:\Bigr{.}$		(3.3c)
		$\displaystyle\max\limits_{i\in I}{\frac{|\beta_{i}-\bar{\beta}_{i}|}{\tilde{\beta}_{i}}}\leq\Delta,$		(3.3d)
		$\displaystyle\hat{\beta}^{min}_{i}\leq\beta_{i}\leq\hat{\beta}^{max}_{i},\ \forall i\in I,$		(3.3e)
		$\displaystyle\max\limits_{i\in I,j\in C_{i}}{\frac{|\gamma_{ij}-\bar{\gamma}_{ij}|}{\tilde{\gamma}_{ij}}}\leq\Delta,$		(3.3f)
		$\displaystyle\hat{\gamma}^{min}_{ij}\leq\gamma_{ij}\leq\hat{\gamma}^{max}_{ij},\ \forall i\in I,j\in C_{i}\Bigl{.}\Bigr{\}},$		(3.3g)

where $\hat{\beta}^{min}_{i}=\min_{s\in\{1,\ldots,S_{i}^{RO}\}}\{\hat{\beta}_{is}\}$, $\hat{\beta}^{max}_{i}=\max_{s\in\{1,\ldots,S_{i}^{RO}\}}\{\hat{\beta}_{is}\}$, $\hat{\gamma}^{min}_{ij}=\min_{s\in\{1,\ldots,T_{ij}^{RO}\}}\{\hat{\gamma}_{ijs}\}$ and $\hat{\gamma}^{max}_{ij}=\max_{s\in\{1,\ldots,T_{ij}^{RO}\}}\{\hat{\gamma}_{ijs}\}$.

Note that the constraints (3.3d) and (3.3f) in the above construction ensure that the maximum normalized deviation of the uncertain parameter $\beta_{i}$ and $\gamma_{ij}$, respectively, from the nominal values among all the products remains under a certain positive scalar $\Delta$. The parameter $\Delta$ is called budget of uncertainty [12]. Constraintd (3.3e) and (3.3g) also define bounds on individual $\beta_{i}$ and $\gamma_{ij}$, respectively.

The budget of uncertainty constraint could also be specified by imposing a limit on the total amount of deviation across all products [13] as follows:

$$\sum\limits_{i\in I}{\frac{|\beta_{i}-\bar{\beta}_{i}|}{\tilde{\beta}_{i}}}+\sum_{i\in I}{\sum_{j\in C_{i}}{\frac{|\gamma_{ij}-\bar{\gamma}_{ij}|}{\tilde{\gamma}_{ij}}}}\leq\Delta.$$

We therefore reformulate the uncertainty set as a polytope:

$$\begin{array}[]{rll}\mathcal{RU}_{\Delta}=\Bigl{\{}\left[\begin{array}[]{c}\boldsymbol{\beta}\\ \boldsymbol{\gamma}\\ \boldsymbol{\sigma}^{\beta}\\ \boldsymbol{\sigma}^{\gamma}\end{array}\right]&\in\mathbb{R}^{2n+2k}:\Bigr{.}&\\ &\sigma_{i}^{\beta}\leq\Delta,&\forall i\in I,\vspace{2mm}\\ &\sigma_{ij}^{\gamma}\leq\Delta,&\forall i\in I,j\in C_{i},\vspace{2mm}\\ &-\tilde{\beta}_{i}\cdot\sigma_{i}^{\beta}\leq\beta_{i}-\bar{\beta}_{i}\leq\tilde{\beta}_{i}\cdot\sigma_{i}^{\beta},&\forall i\in I,\vspace{2mm}\\ &\hat{\beta}^{min}_{i}\leq\beta_{i}\leq\hat{\beta}^{max}_{i},&\forall i\in I,\\ &-\tilde{\gamma}_{ij}\cdot\sigma_{ij}^{\gamma}\leq\gamma_{ij}-\bar{\gamma}_{ij}\leq\tilde{\gamma}_{ij}\cdot\sigma_{ij}^{\gamma},&\forall i\in I,j\in C_{i},\vspace{2mm}\\ &\hat{\gamma}^{min}_{i}\leq\gamma_{ij}\leq\hat{\gamma}^{max}_{ij},&\forall i\in I,j\in C_{i}\Bigl{.}\Bigr{\}}\end{array}$$

(3.4)

Note that the right-hand side of the constraint (3.2a) constitutes an optimization problem. We substitute the uncertainty set $\mathcal{U}$ in  (3.2a) with a hyperplane derived from (3.4), and yield the following optimization problem:

	$\displaystyle\min\limits_{\beta,\gamma,\sigma^{\beta},\sigma^{\gamma}}\quad$	$\displaystyle\sum\limits_{i\in I}{\Bigl{(}-\beta_{i}x_{i}+\sum_{j\in C_{i}}{\gamma_{ij}y_{ij}}\Bigr{)}}$
	$\displaystyle\mathop{\rm s.t.}\quad$	$\displaystyle\sigma_{i}^{\beta}\leq\Delta,$	$\displaystyle\forall i\in I,$		(3.5a)
		$\displaystyle\sigma_{ij}^{\gamma}\leq\Delta,$	$\displaystyle\forall i\in I,j\in C_{i},$		(3.5b)
		$\displaystyle\beta_{i}-\bar{\beta}_{i}\leq\tilde{\beta}_{i}\cdot\sigma_{i}^{\beta},$	$\displaystyle\forall i\in I,\vspace{2mm}$		(3.5c)
		$\displaystyle-\tilde{\beta}_{i}\cdot\sigma_{i}^{\beta}\leq\beta_{i}-\bar{\beta}_{i},$	$\displaystyle\forall i\in I,\vspace{2mm}$		(3.5d)
		$\displaystyle\hat{\beta}^{min}_{i}\leq\beta_{i}$	$\displaystyle\forall i\in I,$		(3.5e)
		$\displaystyle\beta_{i}\leq\hat{\beta}^{max}_{i},$	$\displaystyle\forall i\in I,$		(3.5f)
		$\displaystyle\gamma_{ij}-\bar{\gamma}_{ij}\leq\tilde{\gamma}_{ij}\cdot\sigma_{ij}^{\gamma},$	$\displaystyle\forall i\in I,j\in C_{i},\vspace{2mm}$		(3.5g)
		$\displaystyle-\tilde{\gamma}_{ij}\cdot\sigma_{ij}^{\gamma}\leq\gamma_{ij}-\bar{\gamma}_{ij},$	$\displaystyle\forall i\in I,j\in C_{i},\vspace{2mm}$		(3.5h)
		$\displaystyle\hat{\gamma}^{min}_{i}\leq\gamma_{ij}$	$\displaystyle\forall i\in I,j\in C_{i}$		(3.5i)
		$\displaystyle\gamma_{ij}\leq\hat{\gamma}^{max}_{ij},$	$\displaystyle\forall i\in I,j\in C_{i}$		(3.5j)
		$\displaystyle\sigma_{i}^{\beta},\sigma_{ij}^{\gamma}\geq 0,$	$\displaystyle i\in I,j\in C_{i}.$

By introducing dual multipliers $\mu^{\beta}_{i}$, $\mu^{\gamma}_{ij}$, $\pi_{1i}^{\beta}$, $\pi_{2i}^{\beta}$, $\lambda_{1i}^{\beta}$, $\lambda_{2i}^{\beta}$, $\pi_{1ij}^{\gamma}$, $\pi_{2ij}^{\gamma}$, $\lambda_{1ij}^{\gamma}$ and $\lambda_{2ij}^{\gamma}$ associated with constraints (3.5a)-(3.5j), respectively, the dual problem becomes:

$$\begin{array}[]{rll}\max&\sum\limits_{i\in I}{\bigl{(}-\Delta\mu_{i}^{\beta}-\bar{\beta}_{i}\pi_{1i}^{\beta}+\bar{\beta}_{i}\pi_{2i}^{\beta}+\hat{\beta}^{min}_{i}\lambda_{1i}^{\beta}-\hat{\beta}^{max}_{i}\lambda_{2i}^{\beta}\bigr{)}}&\\ &+\sum\limits_{i\in I}{\sum\limits_{j\in C_{i}}{\Bigl{(}-\Delta\mu_{ij}^{\gamma}-\bar{\gamma}_{i}\pi_{1ij}^{\gamma}+\bar{\gamma}_{i}\pi_{2ij}^{\gamma}+\hat{\gamma}^{min}_{ij}\lambda_{1ij}^{\gamma}\Bigl{.}}}\\ &\hskip 54.06023pt-\hat{\gamma}^{max}_{ij}\lambda_{2ij}^{\gamma}\Bigr{.}\Bigl{)}&\vspace{2mm}\\ \mathop{\rm s.t.}&-\pi_{1i}^{\beta}+\pi_{2i}^{\beta}+\lambda_{1i}^{\beta}-\lambda_{2i}^{\beta}=-x_{i}&i\in I,\vspace{2mm}\\ &-\mu_{i}^{\beta}+\tilde{\beta}_{i}\pi_{1i}^{\beta}+\tilde{\beta}_{i}\pi_{2i}^{\beta}\leq 0,&i\in I,\vspace{2mm}\\ &-\pi_{1ij}^{\gamma}+\pi_{2ij}^{\gamma}+\lambda_{1ij}^{\gamma}-\lambda_{2ij}^{\gamma}=y_{ij}&i\in I,j\in C_{i},\vspace{2mm}\\ &-\mu_{ij}^{\gamma}+\tilde{\gamma}_{ij}\pi_{1ij}^{\gamma}+\tilde{\gamma}_{ij}\pi_{2ij}^{\gamma}\leq 0,&i\in I,j\in C_{i},\vspace{2mm}\\ &\mu_{i}^{\beta},\pi_{1i}^{\beta},\pi_{2i}^{\beta},\lambda_{1i}^{\beta},\lambda_{2i}^{\beta}\geq 0.&i\in I,\vspace{2mm}\\ &\mu_{ij}^{\gamma},\pi_{1ij}^{\gamma},\pi_{2ij}^{\gamma},\lambda_{1ij}^{\gamma},\lambda_{2ij}^{\gamma}\geq 0.&i\in I,j\in C_{i}.\end{array}$$

(3.6)

In accordance with weak duality, any feasible solution to (3.6) provides a lower bound to (3.4). Therefore, we can enforce (3.2a) by adding the following constraint

	$\displaystyle\eta\leq$	$\displaystyle\sum\limits_{i\in I}{\bigl{(}-\Delta\mu_{i}^{\beta}-\bar{\beta}_{i}\pi_{1i}^{\beta}+\bar{\beta}_{i}\pi_{2i}^{\beta}+\hat{\beta}^{min}_{i}\lambda_{1i}^{\beta}-\hat{\beta}^{max}_{i}\lambda_{2i}^{\beta}\bigr{)}+}$
		$\displaystyle\sum\limits_{i\in I}{\sum\limits_{j\in C_{i}}{\Bigl{(}-\Delta\mu_{ij}^{\gamma}-\bar{\gamma}_{i}\pi_{1ij}^{\gamma}+\bar{\gamma}_{i}\pi_{2ij}^{\gamma}+\hat{\beta}^{min}_{i}\lambda_{1ij}^{\gamma}-\hat{\beta}^{max}_{i}\lambda_{2ij}^{\gamma}\Bigr{)}}}$

where $(\mu^{\beta},\pi_{1\cdot}^{\beta},\pi_{2\cdot}^{\beta},\lambda_{1\cdot}^{\beta},\lambda_{2\cdot}^{\beta},\mu^{\gamma},\pi_{1}^{\gamma},\pi_{2}^{\gamma},\lambda_{1}^{\gamma},\lambda_{2}^{\gamma})$ is feasible solution for (3.6). Therefore, the Max-Min problem can be reformulated as:

	$\displaystyle Z_{\Delta}^{*}=\max$	$\displaystyle\hskip 2.84526pt\eta+\sum\limits_{i\in I}{\alpha_{i}p_{i}}$			(3.7a)
	$\displaystyle\mathop{\rm s.t.}$	$\displaystyle\hskip 2.84526pt\eta\leq\sum\limits_{i\in I}{\bigl{(}-\Delta\mu_{i}^{\beta}-\bar{\beta}_{i}\pi_{1i}^{\beta}+\bar{\beta}_{i}\pi_{2i}^{\beta}+\hat{\beta}^{min}_{i}\lambda_{1i}^{\beta}\bigr{.}}$
		$\displaystyle\hskip 14.22636pt-\hat{\beta}^{max}_{i}\lambda_{2i}^{\beta}\bigl{.}\bigr{)}+\sum\limits_{i\in I}{\sum\limits_{j\in C_{i}}{\Bigl{(}-\Delta\mu_{ij}^{\gamma}-\bar{\gamma}_{i}\pi_{1ij}^{\gamma}+\bar{\gamma}_{i}\pi_{2ij}^{\gamma}\Bigr{.}}}$
		$\displaystyle\hskip 14.22636pt\Bigl{.}+\hat{\gamma}^{min}_{ij}\lambda_{1ij}^{\gamma}-\hat{\gamma}^{max}_{ij}\lambda_{2ij}^{\gamma}\Bigr{)}$			(3.7b)
		$\displaystyle-\pi_{1i}^{\beta}+\pi_{2i}^{\beta}+\lambda_{1i}^{\beta}-\lambda_{2i}^{\beta}=-x_{i}$	$\displaystyle i\in I,\vspace{2mm}$		(3.7c)
		$\displaystyle-\mu_{i}^{\beta}+\tilde{\beta}_{i}\pi_{1i}^{\beta}+\tilde{\beta}_{i}\pi_{2i}^{\beta}\leq 0,$	$\displaystyle i\in I,\vspace{2mm}$		(3.7d)
		$\displaystyle-\pi_{1ij}^{\gamma}+\pi_{2ij}^{\gamma}+\lambda_{1ij}^{\gamma}-\lambda_{2ij}^{\gamma}=y_{ij}$	$\displaystyle i\in I,j\in C_{i},\vspace{2mm}$		(3.7e)
		$\displaystyle-\mu_{ij}^{\gamma}+\tilde{\gamma}_{ij}\pi_{1ij}^{\gamma}+\tilde{\gamma}_{ij}\pi_{2ij}^{\gamma}\leq 0,$	$\displaystyle i\in I,j\in C_{i},\vspace{2mm}$		(3.7f)
		$\displaystyle\text{\eqref{prob:reform2:const1} - \eqref{prob:reform2:const4}},$			(3.7g)
		$\displaystyle\mu_{i}^{\beta},\pi_{1i}^{\beta},\pi_{2i}^{\beta},\lambda_{1i}^{\beta},\lambda_{2i}^{\beta}\geq 0.$	$\displaystyle i\in I,\vspace{2mm}$		(3.7h)
		$\displaystyle\mu_{ij}^{\gamma},\pi_{1ij}^{\gamma},\pi_{2ij}^{\gamma},\lambda_{1ij}^{\gamma},\lambda_{2ij}^{\gamma}\geq 0.$	$\displaystyle i\in I,j\in C_{i}.$		(3.7i)

The above formulation constitutes the RO model associated with a budget of uncertainty parameter $\Delta$. This RO model is a conventional linear programming model that can be readily solved using standard linear programming solvers. We report the results of our numerical experiments based on this model in the subsequent section.

We tested the proposed approach on 40 different instances from the data of a major food distributor in the United States. The company sells products to its customers via its geographically scattered divisions. Table 1 lists the name of the test instances (such that each instance represents a specific division), the number of products and number of cross terms (nonzero $\gamma_{ij}$ terms) within a given time period in each instance.

The value of $\Delta$ plays a crucial role in our tests. We select and test the values as follows. First, let us revisit the constraint involving $\Delta$:

	$\displaystyle\max_{i\in I}\left\{\frac{|\beta_{i}-\bar{\beta}_{i}|}{\tilde{\beta}_{i}}\right\}\leq\Delta,\quad\forall i\in I,$
	$\displaystyle\max\limits_{i\in I,j\in C_{i}}{\frac{|\gamma_{ij}-\bar{\gamma}_{ij}|}{\tilde{\gamma}_{ij}}}\leq\Delta,\quad\forall i\in I,j\in C_{i}.$

Obviously, the minimum value for $\Delta$ is zero, in which case all the $\beta_{i}$ values will be equal to their nominal value $\bar{\beta}_{i}$. We can also find the maximum value for $\Delta$ by:

$$\Delta_{max}^{\beta}=\max_{i\in I}\left\{\max\left\{\frac{|\hat{\beta}^{min}_{i}-\bar{\beta}_{i}|}{\tilde{\beta}_{i}},\frac{|\hat{\beta}^{max}_{i}-\bar{\beta}_{i}|}{\tilde{\beta}_{i}}\right\}\right\},$$

and

$$\Delta_{max}^{\gamma}=\max_{i\in I,j\in C_{i}}\left\{\max\left\{\frac{|\hat{\gamma}^{min}_{i}-\bar{\gamma}_{i}|}{\tilde{\gamma}_{i}},\frac{|\hat{\gamma}^{max}_{i}-\bar{\gamma}_{i}|}{\tilde{\gamma}_{i}}\right\}\right\}.$$

Therefore, we choose $\Delta_{max}=\max\{\Delta_{max}^{\beta},\Delta_{max}^{\gamma}\}$. We then divide the interval [0,$\Delta_{max}$] into 10 equally distanced segments, namely $[\Delta_{1},\Delta_{2}]$, $[\Delta_{2},\Delta_{3}],\ldots,[\Delta_{10},\Delta_{11}]$ where $\Delta_{1}=0$ and $\Delta_{11}=\Delta_{max}$, and test each of the $\Delta$ values in our subsequent experiments.

Other parameters of the model are $S_{i}^{RO}$ and $T_{ij}^{RO}$, the numbers of samples used for calculating the means and variances for $\beta_{i}$ and $\gamma_{ij}$. We set these two parameters to 50 and 5, respectively, reflecting the fact that the cross price elasticities are usually associated with smaller intervals compared with the own price elasticities, as found in our real data sets.

To present our results, we show how changing the value of $\Delta$ affects the optimal value of the robust counterpart problem (3.7). To do this, we calculate the percent difference between the maximum and the minimum value of $Z^{*}$ associated with values of $\Delta$ in problem (3.7). The maximum and minimum values of $Z^{*}$ are attained when $\Delta=0$ and $\Delta=\Delta_{max}$, i.e., the least conservative case and the most conservative case, respectively. We summarize the results in Table 2 by only displaying percent of decline from the least conservative to the most conservative case. We observe that the differences are all within in the range $(0.05\%,6.96\%)$ with an average of 1.14%.

The distribution of the percent differences is depicted in the box plot of Figure 2. The detailed results are shown in the Appendix (Table 5).

We further validate our RO approach by comparing results from a set of Monte Carlo simulation tests with our RO solutions. We run a Monte Carlo simulations by taking $S^{sim}$ uniformly distributed samples $\beta_{is}$ from the interval $[\hat{\beta}^{min}_{i},\hat{\beta}^{max}_{i}]$ and $\gamma_{ijs}$ from $[\hat{\gamma}_{ij}^{min},\hat{\gamma}_{ij}^{max}]$, and then solving each of the price optimization problems with given sampling values for $\beta$ and $\gamma$:

$$\begin{array}[]{rl}f(\boldsymbol{\beta}_{s},\boldsymbol{\gamma}_{s})=\max\limits_{p,x,y}&\sum\limits_{i\in I}{(\alpha_{i}p_{i}-\beta_{is}x_{i}+\sum_{j\in C_{i}}{\gamma_{ij}y_{ij}})\vspace{1mm}}\\ \mathop{\rm s.t.}&\eqref{prob:reform2:const1}-\eqref{prob:reform2:const4}.\end{array}$$

(4.1)

The value of the simulation will then be the mean of optimal values from each replication, i.e. $Z_{sim}^{*}=\frac{1}{S^{sim}}\sum\limits_{s}{f(\boldsymbol{\beta}_{s},\boldsymbol{\gamma}_{s})}$.

In these experiments, with $S^{sim}=50$, we compare the values of the simulation $Z^{*}_{sim}$ with the values of RO formulation (3.7). This comparison discloses the degree of conservativeness of the RO solutions. To perform such a comparison, we compare $Z^{*}_{sim}$ with the least conservative case (achieved at $\Delta_{1}$), the most conservative case $\Delta_{max}$ and the average across all $\Delta$ values, and calculate the corresponding difference ratios. For instance the difference ratio in the column under $\Delta_{1}$ is attained from

$$\frac{Z^{*}_{sim}-Z^{*}_{\Delta_{1}}}{Z^{*}_{\Delta_{1}}}.$$

Table 3 summarizes these comparisons. while a negative value of difference ratio in the table implies that the RO value is better than the simulation, a positive value shows conversely. Note that most of the instances (35 out of 40) have positive values in $\Delta_{max}$ column, indicating RO results are worse than the simulation values. However, in most of such cases (30 out of 35) the difference ratios are less than 0.5%, and there are only three cases where the difference ratios are over 1%. This implies that for the most conservative scenarios, RO solution yields results with little losses as compared to the solutions from the simulation. On the other hand, in the least conservative case (under column $\Delta_{1}$), most (31 of 40) differences are negative, demonstrating RO results are better than those from the simulations, while there are only nine instances where it is not and only in one of which the difference ratio is above 0.32%.