Analysis of a chance-constrained new product risk model with multiple customer classes

Accepted Manuscript

Analysis of a chance-constrained new product risk model with multiple customer classes Belleh Fontem, Jeremiah Smith PII: DOI: Reference:

S0377-2217(18)30664-7 https://doi.org/10.1016/j.ejor.2018.07.042 EOR 15280

To appear in:

European Journal of Operational Research

Received date: Revised date: Accepted date:

4 September 2017 7 June 2018 29 July 2018

Please cite this article as: Belleh Fontem, Jeremiah Smith, Analysis of a chance-constrained new product risk model with multiple customer classes, European Journal of Operational Research (2018), doi: https://doi.org/10.1016/j.ejor.2018.07.042

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights • We consider a cost-minimization nonlinear program featuring multiplicative noise • We derive a closed-form expression for the optimal solution in a special case

CR IP T

• For a related special case, we reduce the problem to a sequence of linear programs

• For other cases, we propose a heuristic and compare it to sample average approxi- mation

AC

CE

PT

ED

M

AN US

• Numerical experiments illustrate how various parameters influence the optimal cost

1

ACCEPTED MANUSCRIPT

Analysis of a chance-constrained new product risk model with multiple customer classes Belleh Fontema , Jeremiah Smitha

CR IP T

a College of Business, University of Mary Washington 1301 College Avenue, Fredericksburg, VA 22401, United States. Corresponding Author: Belleh Fontem; E-mail: [email protected]; Phone: +1-540-654-1729; Fax: +1-540-654-2430.

Abstract

AN US

We consider the non-convex problem of minimizing a linear deterministic cost objective subject to a probabilistic requirement on a nonlinear multivariate stochastic expression attaining, or exceeding a given threshold. The stochastic expression represents the output of a noisy system featuring the product of mutually-independent, uniform random parameters each raised to a linear function of one of the decision vector’s constituent variables. We prove a connection to (i) the probability measure on the superposition of a finite collection of uncorrelated expo-

M

nential random variables, and (ii) an entropy-like affine function. Then, we determine special cases for which the optimal solution exists in closed-form, or is accessible via sequential linear

ED

programming. These special cases inspire the design of a gradient-based heuristic procedure that guarantees a feasible solution for instances failing to meet any of the special case conditions. The application motivating our study is a consumer goods firm seeking to cost-effectively man-

PT

age a certain aspect of its new product risk. We test our heuristic on a real problem and compare its overall performance to that of an asymptotically optimal Monte-Carlo-based method called

CE

sample average approximation. Numerical experimentation on synthetic problem instances sheds light on the interplay between the optimal cost and various parameters including the probabilistic requirement and the required threshold.

AC

Keywords: Stochastic programming, Chance-constrained optimization, Multiplicative uncertainty, Risk management

Preprint submitted to European Journal of Operational Research

August 3, 2018

ACCEPTED MANUSCRIPT

1. Introduction 1.1. Background and problem description Traditional mathematical optimization models conveniently ignore the possibility of any un-

CR IP T

certainty in the parameters influencing the optimality or feasibility of the decision vectors involved. In practice, this implicit assumption of perfect knowledge is often a strong enough one to significantly diminish the real-world fidelity of optimal solutions for such classical models.

Chance-constrained optimization has thus emerged as an approach for coping with the nuisance

of parameter uncertainty in optimization problems. The goal of a chance-constrained model is to optimize a function subject to one or more chance constraints which demand explicit probabilis-

AN US

tic guarantees on the satisfaction of linear or nonlinear inequalities featuring noisy parameters. That is, chance-constrained models aspire not just to optimize performance, but also to be robust in the face of uncertainty. If, in the absence of uncertainty, either the objective function,

or at least one member of the set of constraint functions is nonlinear, then we have a nonlinear program, and by infecting any of its constraint parameters with randomness, we turn the optimization problem into a chance-constrained nonlinear program (CCNLP).

M

The discipline of chance-constrained optimization traces its origins to the pioneering papers of Charnes et al. (1958), and Charnes and Cooper (1959) both of which addressed stochastic

ED

problems in heating oil inventory planning. Over the years, the body of theoretical and computational knowledge has grown severalfold with notable contributions being those of Pr´ekopa (1973, 1993, 1995); Pr´ekopa et al. (1998, 2011), Sepp¨al¨a (1972), Dentcheva et al. (2000), Nemirovski

PT

and Shapiro (2006, 2007), Luedtke et al. (2010) to name a few (see Geletu et al. (2013) for a more comprehensive survey). More recently, techniques developed to solve chance-constrained

CE

problems have become increasingly popular among engineers and decision-makers as tools for handling uncertainty in various application domains such as power systems management (Kargarian et al., 2016), robotic path control (Blackmore et al., 2011), network machine learning

AC

(Doppa et al., 2010), portfolio optimization (Fakhar et al., 2018), and forest disaster response (Wei et al., 2015). However, coming up with exact algorithms for chance-constrained models is quite a chal-

lenging task principally because (i) chance constraints are in general non-convex, and (ii) it is typically not the case that we have a closed-form probability formula at hand to a priori check the feasibility of any given candidate solution (Pr´ekopa, 1995). Hence, the absence of generaliz3

ACCEPTED MANUSCRIPT

able exact methods frequently compels us to resort to heuristic alternatives such as back-mapping (Wendt, et al., 2002; Li et al., 2008; Zhang and Li, 2011), or to asymptotically optimal MonteCarlo-based methods like sample average approximation (Luedtke and Ahmed, 2008; Pagnoncelli et al., 2009). Nevertheless, special chance-constrained linear programs (such as those with

CR IP T

Gaussian uncertainty and known covariances) do possess the desirable properties of chance con-

straint convexity, and easy a priori feasibility checking (Ahmed and Shapiro, 2008; Zhang et al., 2016). An interesting question therefore, is whether it is possible for a CCNLP to be endowed

with at least one of the aforementioned desirable properties. The answer to that question might be useful in exploring for instance, a (generally non-convex) CCNLP of the form:

x

AN US

X   z x? , min cn xn

(1)

n∈N

subject to     Y   − ak+n xn  (e an ) n ≥ b ≥ p Pr 1 −

(2)

n∈N

(3)

M

xn ≥ ln for all n ∈ N

ED

where b ∈ (0, 1), p ∈ (0, 1), cn > 0, kn ∈ (0, 1], a+n > 1, ln ≥ 0 are deterministic scalars for all X n ∈ N and |N| ≥ 1, kn = 1. The constrained probability measure (2) is our chance constraint,

PT

n∈N  and the quantity e an ∈ A , e a1 ,e a2 , . . . ,e a|N| is a continuous random variable that is uniformly  distributed on the support 0, a+n . To keep the analysis tractable, we assume mutual indepen        dence of all elements in A, thereby implying that Cov e an ,e am , E e ane am − E e an E e am = 0

for all n , m, and we justify the imposition of a continuous uniform distribution on e an ∈ A

CE

because of the continuous uniform distribution’s maximum differential entropy property. The

e (distributed on continuous domain Ω) with density differential entropy of a random variable ω R function fωe (·) is defined as H (e ω) , − Ω fωe (ω) log ωdω, and roughly speaking, inversely re-

AC

e. Being the distribution lates to the amount of prior information that is intrinsically built into ω

with the greatest differential entropy, the continuous uniform distribution therefore assumes the

least amount of prior information when serving as a model for the uncertainty of any unknown continuous quantity with known bounds.

4

ACCEPTED MANUSCRIPT

1.2. Motivating application and research contributions The motivation for studying Problem (1) arises out of a certain consumer goods firm’s new   Y − ak+n xn (e product risk model in which the expression 1 − an ) n ∈ [0, 1) represents a fractional and n∈N

stochastic index measure reflecting the market’s aggregate perception of the attractiveness (i.e.,

CR IP T

desirability) of a new product in the firm’s portfolio. The product is marketed at |N| disparate

customer classes, and its aggregate desirability in the market (measured via survey response data)

is taken to depend on two sets of ingredients. The first set is stored in a decision vector of real  0 numbers x = x1 x2 . . . x|N| in the non-negative orthant R|N| + for which xn ≥ 0 represents

a firm-provided incentive that entices class n ∈ N members to promote the new product on

AN US

the firm’s behalf. Furthermore, the unit cost of an incentive destined for class n members is

cn > 0 and out of strategic or operational necessity, the firm is constrained by a lower bound  ln ≥ 0 on each decision variable xn . The second set of ingredients is a collection kn , a+n ,e an n∈N of deterministic and stochastic parameters which govern how members of the various customer classes typically react to the firm’s incentives.

More concretely, the decision variable xn ≥ 0 models the time and effort expended by the

M

firm to induce class n ∈ N customers to effectively act as a surrogate arm of the firm’s salesforce

whereas kn , a+n and e an dictate the predictable and random aspects relating to how sensitive class  Y − ak+n xn (e an ) n therefore maps an input vector x of costly n customers are to xn . The model 1 −

ED

n∈N

firm-provided promotion incentives to a random (and hence risky) output fraction that measures the new product’s overall desirability in the market.

PT

Expressions such as the posited stochastic model 1 −

Y

n∈N

(e an )

  − ak+n xn n

are effectively probabilis-

tic versions of models belonging to the family of so-called production functions extensively used

CE

in economics (Shahbaz et al., 2017; Grassetti and Hunanyan, 2018), agronomy (Araya et al., 2018; Zhang et al., 2018) and engineering (Wibe, 1984). At heart, a production function1 describes the relationship between a valuable output and a set of controllable input factors which

AC

collectively determine the output quantity (e.g., an agricultural production function could model how a certain crop yield depends on necessary inputs like labor, capital and technology). To capture economic reality, a production function must have the attributes of essentiality,

monotonicity, upper semi-continuity, attainability, and quasi-concavity in every input. In what 1 For

a detailed text on production functions, we refer the reader to Bairam (1994).

5

ACCEPTED MANUSCRIPT

immediately follows, we show that the conditional value of 1 −

Y

(e an )

n∈N

  − ak+n xn n

qualifies as a

production function.  0 Let a , a1 a2 . . . a|N| be a vector of hypothesized values of the respective random  parameters in A = e a1 ,e a2 , . . . ,e a|N| where 0 < an ≤ a+n , an , 1 for all n ∈ N. Now, suppose that n∈N

• Essentiality: PDes vanishes in the absence of any inputs i.e., −

(an )

n∈N



kn + an

 (0)

= 1 − 1|N| = 1 − 1 = 0.

AN US

PDes (0|a) = 1 −

Y

CR IP T

0 ∈ R|N| is a vector of zeros and define the mapping PDes : R|N| + 7−→ [0, 1) to be PDes (x|a) ,  Y − ak+n xn (an ) n . Then, it is clear that PDes has the following properties. 1−

• Strict monotonicity (and therefore, monotonicity): PDes is strictly increasing in all inputs i.e.,

M

    Y  ∂PDes (x|a) ∂  − ak+n xn  1 − (an ) n  = ∂xn ∂xn n∈N     Y kn − ak+n xn − akm+ xm n (am ) m = + (an ) log (an ) > 0 for all xn ≥ 0, n ∈ N. an m∈N

ED

m,n

PT

• Continuity (and therefore, upper semi-continuity): PDes is smooth all over its domain i.e.,

CE

t     !t Y ∂ PDes (x|a) − kn x − km x t−1 kn (an ) a+n n logt (an ) (am ) a+m m > 0 for all xn ≥ 0, n ∈ N and t = 1, 2, . . . = (−1) t + ∂xn an m∈N m,n

AC

• Attainability: it is always possible meet any value on PDes ’s range by proportionally and sufficiently increasing all inputs i.e.,

6

ACCEPTED MANUSCRIPT

PDes (x|a) = 1 −

Y

−

(an )

n∈N



kn a+ n



xn

> 0 =⇒ lim PDes (λx|a) = lim 1 − λ↑∞

where λ > 0 and sup PDes (x|a) = lim 1 − xn ↑∞

x∈R|N| +

λ↑∞

Y

−

(an )

n∈N



kn a+ n



xn

Y

−

(an )

n∈N



kn a+ n

 λxn

= sup PDes (x|a) = 1 x∈R|N| +

= 1 for all n ∈ N.

CR IP T

                

• Strict concavity (and therefore quasi-concavity) in every input: PDes obeys the law of diminishing marginal returns with respect to every input i.e.,

AN US

    Y  ∂2 PDes (x|a) ∂2  − ak+n xn  (an ) n  = 2 1 − 2 ∂xn ∂xn n∈N     !2 Y kn − ak+n xn − akm+ xm 2 n (am ) m log (an ) < 0 for all xn ≥ 0, n ∈ N. = − + (an ) an m∈N m,n

eDes (x) , 1− Hence, for a given new product in the firm’s portfolio, the model P

Y

n∈N

(e an )

  − ak+n xn n

M

is a stochastic production function2 that translates the firm’s incentivization strategy x into a frac-

ED

tional desirability measure for that product. In the firm’s context, the uniform random variable  e an ∈ 0, a+n (where a+n > 1) is called the propensity of class n because the magnitude of its logarithm represents the uncertain (i.e., noisy) component of the degree to which a typical class n member will react to infinitesimal changes in the firm’s incentives   !(when these incentives are still

PT

an ) low). To see how this is the case, note that lim Pr (e xn ↑∞

−

kn a+ n

xn

= 0 = 1 almost surely, but for low

  − ak+n xn

  ≈ 1 − ak+nn xn log (e an ), and taking   k the first derivative with respect to xn yields the gradient magnitude a+nn log (e an ) . The random   an ) . In contrast, the dimensionvariable log (e an ) is then the stochastic component of ak+nn log (e n

CE

an ) values of xn , we obtain the first order approximation (e

kn a+n

> 0 (with kn being a relative weight) is called the responsivity of class n, and is   the deterministic component of the gradient magnitude ak+nn log (e an ) . Essentially, in the firm’s

AC

less ratio

2 Insights from the analysis of P eDes are potentially useful in broader contexts involving the investigation of similar eDes might bear on stochastic versions production functions under uncertainty. For example, results from our treatment of P Q
of Johann von Th¨unen’s exponential production function P JVT = A 3n=1 1 − exp (−ˆan φn ) which closely resembles PDes (x|a). In von Th¨unen’s production function, P JVT is some agricultural output with supremum A > 0, the variables φ1 ≥ 0, φ2 ≥ 0, φ3 ≥ 0 are the respective inputs of labor, capital, and fertilizer whereas aˆ 1 > 0, aˆ 2 > 0, aˆ 3 > 0 are instantaneous rates of decline of the marginal productivity of labor, capital, and fertilizer respectively (Humphrey, 1997; Mishra, 2010).

7

ACCEPTED MANUSCRIPT

multiplicative sensitivity model,

kn a+n

approximates the predictable element pertaining to how sen-

sitive class n customers are to its incentives, whereas log (e an ) is an approximate description of the probabilistic component of that sensitivity.

CR IP T

At this juncture, it is important to mention that the point of discussing the first-order ap    − ak+n xn an ) n proximation (e ≈ 1 − ak+nn xn log (e an ) in the preceding paragraph was to introduce and

explain the genesis of the terms “propensity” and “responsivity” rather than to suggest that a

first-order approximation is sufficiently powerful to capture all the model’s moving parts (in fact,   Y − ak+n xn (e 1− an ) n involves an infinite power series, and so employing a first-order approximation n∈N

of each term (e an )

  − ak+n xn

would coarsen the model). Note that a class’ maximum propensity value p   a+ an = 2 √n 3 ) inversely relates to the class’ responsivity a+n (which is linearly proportional to Var e

AN US

kn a+n .

n

This inverse relationship models the firm’s empirical data which suggests that greater un-

certainty over the true value of e an is plausibly explainable in large part by weaker deterministic

sensitivity to changes in the incentive level.

The firm’s objective−for a specific new product in its portfolio−is to find the least costly

feasible incentivization strategy x? capable of ensuring that the product in question will enjoy

M

an aggregate desirability index of at least b ∈ (0, 1) with minimum reliability (i.e., probability)

p ∈ (0, 1). Equivalently, the firm’s goal is to select x? such that the risk of violating (2) is at most

ED

1 − p where p is a chosen probabilistic feasibility standard. To ease the notation, we introduce

for class n’s responsivity, the symbol

kn =⇒ rn ∈ (0, 1) since kn ∈ (0, 1], a+n > 1 for all n ∈ N, and a+n rn rn 1 = + =⇒ ∈ (0, 1) since a+n > 1 for all n ∈ N. kn an kn

CE

PT

rn ,

We confirm that (2) is well-behaved because for all n ∈ N, 0 ≤ xn < ∞        −rn xn Y Y      (e ) (e ) r a log a  n n n −rm xm −rn xn  1 −   > 0   (e ) (e ) a > 0 a > 0 = Pr    n m      cn m∈N n∈N m,n   Y   ∂ −rn xn  (e Pr 1 − =⇒ an ) > 0 > 0 ∂cn xn n∈N

AC

  ∂ + an > 1 =⇒ Pr  ∂cn xn

8

ACCEPTED MANUSCRIPT

and   Y   −rn xn  (e an ) lim cn xn = 0 ⇐⇒ lim Pr 1 − = 0 = 1 x1 ↓0,...,x|N| ↓0 x1 ↓0,...,x|N| ↓0 n∈N   n∈N Y X   (e an )−rn xn = 1 = 1 a.s. lim cn xn = ∞ ⇐⇒ lim Pr 1 − x1 ↑∞,...,x|N| ↑∞ x1 ↑∞,...,x|N| ↑∞ X

n∈N

n∈N

CR IP T

                  

Verbally, the probability of an arbitrary incentivization strategy x being feasible for (2) is monotonically increasing with respect to each cost component of that strategy. Also, a zero cost strategy will always violate (2) whereas an infinitely costly one will satisfy (2) almost surely. In

this paper, Proposition 5 provides a closed-form expression for the left-hand side of (2), thereby

AN US

showing that Problem (1) belongs to the family of chance-constrained programs with easy a

priori feasibility checking. We thus give a positive answer to the question raised in the third paragraph of Subsection 1.1. More broadly, we offer the following four contributions to the literature:

• Given a vector of decisions x, Proposition 1 asserts that feasibility for Problem (1) is

M

strongly connected to the probability that the sum of mutually-uncorrelated exponential random variables (collectively dependent on x) does not exceed a certain entropy-like

ED

affine function evaluation of x.

• Under two restrictive assumptions, we derive a closed-form expression for Problem (1)’s optimal solution (Proposition 2), and for the chance constraint’s shadow price (Corollary

PT

1). We then weaken one of the assumptions to usher in a regime in which Problem (1) is amenable to a linear programming solution framework (Proposition 3). We show that

CE

an important parameter in that framework is the lone positive root of a certain quadratic equation (Proposition 4). To summarize, we identify special (albeit restrictive) cases in which we can quickly and exactly determine Problem (1)’s optimal solution. These special

AC

cases are valuable for two reasons. First, they show that at least for the class of problems of the form of Problem (1), there are situations in which we could do better than resorting to (asymptotically optimal) Monte-Carlo-based solution techniques. Second, they provide useful information for designing heuristic procedures aimed at pragmatically addressing Problem (1) in very general situations.

• Building on the preceding contribution, we exploit the result in Corollary 1 to design 9

ACCEPTED MANUSCRIPT

a heuristic that covers cases outside the jurisdiction of the previously-mentioned special cases. The proposed heuristic begins by solving a special case problem that yields a “warm start” initial infeasible solution. Then, using estimated shadow price information at each iteration, the heuristic gradually works its way toward the feasible space with the goal of

CR IP T

minimizing the journey’s cumulative cost. Using real parameter settings from a consumer

goods firm, we compare our heuristic procedure to the Monte-Carlo-based sample average approximation approach of Luedtke and Ahmed (2008).

• We artificially generate problem instances that conform to the special cases. Then, we experiment on these instances so as to validate the theoretical analysis and to better under-

AN US

stand how the optimal cost reacts to changes in the propensity regime.

Outline. The organization of this paper is as follows. In Section 2, we analyze the problem theoretically and derive some structural properties. Section 3 then provides a special case solution algorithm that exploits these properties. In Section 4, we propose a heuristic procedure to cover instances beyond our special case algorithm’s reach, then we derive a mixed integer lin-

M

ear programming problem that incarnates an asymptotically optimal solution technique known as sample average approximation. Section 5 experimentally verifies the findings on a real problem instance, as well as on numerous replications of randomly-generated problem instances. In

ED

Section 6, we round up the discussion and offer possible fruitful extensions of our work. Notation. Throughout this paper, lowercase boldface shall denote vectors, the tilde (∼) su-

PT

perscript over a letter will indicate randomness, and on its own, will be an abbreviation for the   phrase “is distributed as”. For instance, e θn (xn ) ∼ Exp rn1xn means that the random variable

e θn (xn ) is exponentially-distributed with rate parameter

1 rn x n .

We will use the lowercase variable t

CE

either as a local operating index, or as the variable of a Riemann integral.

AC

2. Main results

We begin by characterizing feasibility for (2) in terms of the probability measure belonging to

the superposition of a set of mutually-uncorrelated exponential random variables. To set the stage for this characterization, we associate to each incentivization strategy x, an index set N>0 (x), a set T (x) of mutually-uncorrelated exponential random variables, and an entropy-like real-valued affine function h(x). Their definitions follow.

10

ACCEPTED MANUSCRIPT

Definition 1. (4) (5)

CR IP T

N>0 (x) , {n ∈ N : xn > 0} ⊆ N n o T (x) , e θn (xn ) : n ∈ N>0 (x) ! h i 1 where e θn (xn ) ∼ Exp for all n ∈ N>0 (x) and Cov e θn (xn ) , e θm (xm ) = rn xn h i h i h i E e θn (xn ) e θm (xm ) − E e θn (xn ) E e θm (xm ) = 0 for all n , m ! X rn h (x) , log(1 − b) − rn xn log . kn n∈N (x) >0

X

(6)

AN US

rn We remark that the affine function in (6) is equivalently h (x) = log(1 − b) − rn xn log k n n∈N

!

rn where x ∈ R|N| + . To boot, b ∈ (0, 1) =⇒ log(1−b) < 0 and for all n ∈ N, we obtain kn ∈ (0, 1) =⇒ ! ! X X rn rn − rn xn log > 0. Therefore, the sign of h (x) = log(1 − b) − rn xn log could kn kn n∈N>0 (x) n∈N>0 (x) be positive, zero, or negative depending on the composition of x. Proposition 1 below establishes

a relationship between (2), superposed mutually-uncorrelated exponential variates, and h(x).

M

Proposition 1 (Chance constraint feasibility characterization) Given a decision vector x, let

ED

e θn (xn ) ∈ T (x) for every n ∈ N>0 (x) be an exponentially-distributed random variable with rate    X  e θn (xn ) ≤ · the cumulative parameter rn1xn . Moreover, denote by FPn∈N (x) eθn (xn ) (·) , Pr  >0 n∈N>0 (x) X e distribution function of θn (xn ). Then, x satisfies (2) if and only if FPn∈N (x) eθn (xn ) (h (x)) is >0

n∈N>0 (x)

PT

at least p. In other words,

CE

  Y   −rn xn  (e an ) ≥ b ≥ p ⇐⇒ FPn∈N (x) eθn (xn ) (h (x)) ≥ p. Pr 1 − >0

(7)

n∈N

Proof. See Appendix A.

AC

Essentially, (7) tells us how to switch between multiplicative and additive uncertainty in the

chance constraint. To see Proposition 1 in action, consider the following toy Example 1. Example 1 Suppose N = {1, 2, 3}, k1 = k2 = k3 = 13 , a+1 = 2, a+2 = 7, a+3 = 5, and let us adopt  the convention that e an ∼ unif 0, a+n means e an ∈ A is uniformly-distributed on the continuous 11

ACCEPTED MANUSCRIPT

 interval 0, a+n where n ∈ N and a+n > 1. Now, define 1

r1 ,

k1 = 3 = 0.167 a+1 2

r2 ,

k2 = 3 = 0.048 a+2 7

r3 ,

k3 = 3 = 0.067 a+3 5

CR IP T

1

1

 and given the column decision vector x , x1

mutually-uncorrelated exponential random variables

x2

0 x3 , define the following family of

AN US

! ! 1 1 = Exp r1 x1 0.167x1 ! ! 1 1 e = Exp θ2 (x2 ) ∼ Exp r2 x2 0.048x2 ! ! 1 1 e θ3 (x3 ) ∼ Exp = Exp r3 x3 0.067x3 e θ1 (x1 ) ∼ Exp

M

where x1 > 0, x2 > 0, x3 > 0. Moreover, for b ∈ (0, 1), define the affine function

ED

         0.167   0.048   0.067   + 0.048x2 log   + 0.067x3 log   . h (x) , log(1 − b) − 0.167x1 log  1 3

1 3

1 3

PT

 i Proposition 1 asserts that the mutual independence of e a1 ∼ unif 0, a+1 = unif (0, 2] ,e a2 ∼  i  i + + unif 0, a2 = unif (0, 7], and e a3 ∼ unif 0, a3 = unif (0, 5] translates to

CE

    a1 )−0.167x1 · (e a2 )−0.048x2 · (e a3 )−0.067x3 ≥ b = Pr e Pr 1 − (e θ1 (x1 ) + e θ2 (x2 ) + e θ3 (x3 ) ≤ h (x)

AC

, Feθ1 (x1 )+eθ2 (x2 )+eθ3 (x3 ) (h (x)) .

12

ACCEPTED MANUSCRIPT

The equivalence relationship in (7) means that we can reformulate Problem (1) to (3) as X   cn xn z x? = min x

(8)

n∈N

FPn∈N

e

θ (x ) >0 (x) n n

CR IP T

subject to −1 (h (x)) ≥ p ⇐⇒ h (x) ≥ FP

xn ≥ ln for all n ∈ N.

−1 where FP

e

n∈N>0 (x) θn (xn )

 (p) , inf h (x) : FPn∈N

tribution (i.e., quantile) function of FPn∈N e

n∈N>0 (x) θn (xn )

e θ (x ) >0 (x) n n

(p)

(9)

(10)

 (h (x)) ≥ p is the inverse cumulative dis-

(·) evaluated at p. We do not have a conX e (p) except when |N>0 (x)| = 1 in which case θn (xn ) = e

θ (x ) >0 (x) n n

AN US

−1 cise expression for FP

e

n∈N>0 (x) θn (xn )

e θn (xn ) is (by definition) exponentially-distributed with rate parameter

1 rn xn

n∈N>0 (x)

and quantile function

Fe−1 (p) = −rn xn log(1 − p) for n ∈ N>0 (x). The function Fe−1 (p) turns out to be quite useful θn (xn )

θn (xn )

because if p lies strictly below a certain critical value, and if a particular subset of {ln }n∈N is

composed of only zero values, then we can use Fe−1 (p) to compute a closed-form solution for θn (xn )

under which the solution is valid.

M

Problem (8). Proposition 2 is a statement of this closed-form solution, and the circumstances

ED

Proposition 2 (Special case closed-form solution) Let crnn be the cost-responsivity ratio for class ( ) cn n ∈ N, and define n , arg min . Now, suppose that n∈N rn

CE

PT

            

p<1−

rn kn

(11)

ln = 0 for all n , n, where n ∈ N.

AC

Then, the components of the optimal solution x? for Problem (8) are            log(1−b)     > 0 , max l   n rn    1    rn log kn · 1−p ? xn =       0

if n = n

(12)

otherwise.

Proof. See Appendix B.

That is, if the minimum required reliability p is small enough, and it is feasible to provide 13

ACCEPTED MANUSCRIPT

zero incentives to any class thatis not n ∈ N, then the firm need only provide class n with an      log(1−b)    . A corollary to Proposition 2 is that under (11), we can express l , incentive of max  n   1  rn log krnn · 1−p  the chance constraint’s shadow price in closed-form. Corollary 1 is a formal statement of this

CR IP T

derivative result.

Corollary 1. Assuming that (11) holds, the optimal differential cost of ensuring probabilistic feasibility (i.e., the chance constraint’s shadow price) is ?


if

log(1−b)   rn 1 rn log kn · 1−p

> ln

(13)

otherwise.

AN US

∂z x ∂p

  cn log(1−b)     >0  −   1  (1−p)3 rn log2 krnn · 1−p =      0

Proof. When (11) holds, (12) is also true, and multiplying both sides of (12) by cn we obtain               (1 c log − b) n   ? ?   c z x = cn xn = max  l , > 0.  n n     rn 1     r log ·  n kn 1−p 

M


(13) is the straightforward result of taking the first derivative of z x? with respect to p.

ED

In Section 4, we will use this shadow price information to help steer a heuristic procedure designed to handle Problem (8) in situations outside the scope of the restrictive cases. At this juncture, one might ask whether it is still possible to find an optimal solution if we generalize

PT

(11) by resetting its second condition to the default ln ≥ 0 for all n ∈ N. Proposition 3 settles this question affirmatively by confirming that we can in fact use linear programming as a solution

CE

method.

AC

  Proposition 3 (Special case linear programming certificate of optimality) Define Λ b, n , ! X rn rn ln log and suppose log(1 − b) − kn n∈N\n p<1−

14

rn . kn

(14)

ACCEPTED MANUSCRIPT

  Then, there exists a constant v? such that Λ b, n ≤ v? < ∞, and the linear program X   cn xn z x? = min x

(15)

n∈N

h (x) ≥ v?

CR IP T

subject to (16)

xn ≥ ln for all n ∈ N has the same optimal solution as that of Problem (8).

AN US

Proof. See Appendix C.

(17)

Remark 1 We observe that (11) is indeed a special case of (14) because when we tighten (17)   by forcing ln = 0 for all n , n (where n ∈ N), the function Λ b, n reduces to log(1 − b), and we

recover (12) from (14) and Problem (15) in the following manner: By (14): p < 1 −

rn kn

ED

M

! ! rn rn rn 1 1 =⇒ ∃xn ∈ [0, ∞) : log(1 − b) ≥ rn xn log · since p < 1 − ⇐⇒ log · <0 kn 1 − p kn kn 1 − p log (1 − b)  . =⇒ ∃xn ∈ [0, ∞) : xn ≥ (18) r 1 rn log knn · 1−p

PT

(17) gives us xn ≥ ln , and combining (18) with the fact that Problem (15) is a minimization

AC

CE

problem, we obtain the familiar closed-form solution

           log(1−b)     > 0 max l ,   n rn   1   rn log kn · 1−p   ? xn =       0

if n = n otherwise.

Proposition 3 provides the assurance that as long as the minimum required reliability p is

sufficiently low, we can solve nonlinear, non-convex Problem (8) by solving instead a linear pro  gram (i.e., Problem (15)) featuring a certain quantity v? whose lower bound is Λ b, n . The exact

determination of v? will certainly be facilitated by knowledge of its upper bound. Proposition 4 below gives us a simple recipe for calculating an upper bound v?U B on v? . 15

ACCEPTED MANUSCRIPT

Proposition 4 (Upper bound on v? ) Suppose p < 1 − u (ˆx) , min x

X

rn kn ,

and let xˆ solve the linear program (19)

cn xn

n∈N

CR IP T

subject to X h (x) ≥ 2 rn ln

(20)

n∈N

xn ≥ ln +

log (1 − b)   for all n ∈ N. r 1 rn log knn · 1−p

The following inequality holds:

v ≤

1−

where

1 αˆ

X

AN US

?



1

(1 − p) αˆ

1+

M

αˆ , 1 +

1 1−p

, v?U B

1 >2 1− p

(22)

(23)

= 0.

ED

is the lone positive root of αˆ 2 − 2αˆ −

rn xˆn

n∈N

s

(21)

Proof. See Appendix D.

The last two propositions show that our CCNLP is conditionally substitutable by a determin-

PT

istic linear program with the same number of constraints 1 + |N| as the CCNLP. To buttress this point, let us explore Example 2 below.

CE

Example 2 Suppose N = {1, 2, 3}, k1 = k2 = k3 = 31 , a+1 = 2, a+2 = 7, a+3 = 5 and c1 = 85, c2 =

AC

60, c3 = 74. By definition, r1 = (

= 0.167, r2 =

k2 a+2

(

= 0.048, r3 =

k3 a+3

= 0.067 and

) 85 60 74 = min , , = min {510, 1260, 1110} = 510 0.167 0.048 0.067 ( ) cn =⇒ n = arg min =1 n∈N rn rn r1 0.167 1 =⇒ 1 − =1− =1− 1 = kn k1 2 3

cn min n∈N rn

)

k1 a+1

16

ACCEPTED MANUSCRIPT

Additionally, let l1 = 29, l2 = 10, l3 = 8 and b = 0.95. Then, ! rn Λ b, n = log(1 − b) − rn ln log kn n∈N\n       0.067    0.048     = −1.211   = log(1 − 0.95) − 0.048(10) log   + 0.067(8) log  X



X

!

1 3

1 3

CR IP T



rn k n n∈N>0 (x)         0.048   0.067    0.167        = log(1 − 0.95) − 0.167x1 log  1  + 0.048x2 log  1  + 0.067x3 log  1  .

h (x) = log(1 − b) −

rn xn log

3

rn kn

=

1 2

3

(say p = 0.45) and h (x) as specified above, note that the

linear program

AN US

For some arbitrary p < 1 −

3

u (ˆx) = min 85x1 + 60x2 + 74x3 x

subject to

h (x) ≥ 2 (0.167(29) + 0.048(10) + 0.067(8)) log (1 − 0.95)   1 0.167 log 0.167 · 1 1−0.45

M

x1 ≥ 29 +

3

PT

ED

log (1 − 0.95)   x2 ≥ 10 + 1 0.167 log 0.167 · 1 1−0.45 3

log (1 − 0.95)   x3 ≥ 8 + 1 0.167 log 0.167 · 1 1−0.45 3

CE

 0  0 yields the minimizer xˆ = xˆ1 xˆ2 xˆ3 = 217.588 198.588 196.588 . Next, introduce q q 1 1 the constant αˆ = 1 + 1 + 1−p = 1 + 1 + 1−0.45 = 2.679 and compute 1−

1 αˆ

X

AC



rn xˆn

n∈N

1

(1 − p) αˆ

= v?U B =



1−

1 2.679



(0.167(217.588) + 0.048(198.588) + 0.067(196.588)) 1

(1 − 0.45) 2.679

= 46.085.

Alongside each other, Propositions 3 and 4 reveal that owing to p = 0.45 < 1 − 17

rn kn

=

1 2,

ACCEPTED MANUSCRIPT

i  i h   a1 ∼ unif 0, a+1 = there exists some v? ∈ Λ b, n , v?U B =⇒ v? ∈ [−1.211, 46.085] such that if e  i  i unif (0, 2] ,e a2 ∼ unif 0, a+2 = unif (0, 7], and e a3 ∼ unif 0, a+3 = unif (0, 5] are mutuallyindependent, then the CCNLP

  z x? = min 85x1 + 60x2 + 74x3

CR IP T

x

subject to   Pr 1 − (e a1 )−0.167x1 · (e a2 )−0.048x2 · (e a3 )−0.067x3 ≥ 0.95 ≥ 0.45 x1 ≥ 29, x2 ≥ 10, x3 ≥ 8

AN US

has an identical optimal solution x? to that of the deterministic linear program   z x? = min 85x1 + 60x2 + 74x3 x

subject to

        0.048   0.067    0.167   + 0.048x2 log   + 0.067x3 log   ≥ v? h (x) = log(1 − 0.95) − 0.167x1 log  1 3

1 3

M

x1 ≥ 29, x2 ≥ 10, x3 ≥ 8

1 3

ED

which happens to possess the same number of constraints as the CCNLP. Propositions 3 and 4 effectively provide most of what we would need to design a solution algorithm for Problem (8) when p < 1 −

rn kn .

Crucially, what they lack is a means of verifying

PT

the feasibility of any arbitrary decision vector x for the chance constraint. Recall that given any

CE

choice of x, Proposition 1 equates the chance constraint with the requirement that the probability X e measure on the sum θn (xn ) of uncorrelated exponential variates evaluated at the output of n∈N>0 (x)

some affine function h (x), be at least p. What we therefore require is a formula for computing

AC

that probability measure FPn∈N

the following.

e

θ (x ) >0 (x) n n

(h (x)) given our selection of x. First, let us formally define

18

ACCEPTED MANUSCRIPT

Definition 2. M (x) , {rn xn : n ∈ N>0 (x)}

(24)

i(n) , min {t : rn xn = rt xt } for all n ∈ N>0 (x)

(25)

I (x) , {i(n) : n ∈ N>0 (x) , i(n) unique}     
 1 if ri(n) xi(n) = rt xt δ ri(n) xi(n) , rt xt ,      0 otherwise

(26)

CR IP T

t∈N>0 (x)

(27)


Mi(n) (x) , rt xt ∈ M (x) : δ ri(n) xi(n) , rt xt = 1 for all i(n) ∈ I (x) .

M (x) =

[

i(n)∈I(x)

and Mi(n) (x)

Mi(n) (x)

\

AN US

Remark 2 Note that for all n, m ∈ N>0 (x),

(28)

Mi(m) (x) = ∅ for all i(n), i(m) ∈ I (x) =⇒ i(n) , i(m).

M

 Essentially, if n ∈ N>0 (x) and i(n) ∈ I (x), then the family Mi(n) (x) i(n)∈I(x) is a disjoint col-

lection of smaller sets in which that labeled Mi(n) (x) is a unique and homogenous set containing

ED

every instance of ri(n) xi(n) in the parent set M (x). That is, for n ∈ N>0 (x) and i(n) ∈ I (x), the integer Mi(n) (x) ≥ 1 represents the multiplicity of element ri(n) xi(n) in M (x). The multiplicity

AC

CE

follows.

PT

of each element in M (x) is central to the computation of FPn∈N

19

e

θ (x ) >0 (x) n n

(h (x)). Proposition 5

ACCEPTED MANUSCRIPT

Proposition 5 (Closed-form expression for the chance constraint’s left-hand side) For n ∈ n o N>0 (x) , i(n) ∈ I (x), and j ∈ 1, . . . , Mi(n) (x) let  ∆i(n), j ,  Mi(n) (x) − j !

lim

s→− r

1 i(n) xi(n)

d|Mi(n) (x)|− j ds|Mi(n) (x)|− j

Then,

  !|Mi(n) (x)| Y  1 1   s +  . (29)  ri(n) xi(n) s + rn1xn  n∈N>0 (x)      Γ j, ri(n)h(x) xi(n)   r x
j 1 − i(n) i(n)  ( j − 1)! 

if h (x) ≥ 0

otherwise.

AN US

 i(n) (x)|  X |MX  ∆i(n), j    Y     i(n)∈I(x) j=1 rn xn FPn∈N (x) eθn (xn ) (h (x)) =    >0 n∈N>0 (x)       0

CR IP T

1

(30)

  R∞ where Γ j, ri(n)h(x) t j−1 exp(−t)dt is the upper incomplete Gamma function. h(x) xi(n) , ri(n) xi(n)

Proof. See Appendix E.

The work in this section has enabled us to discern two special cases for which Problem (8)

M

either has a closed-form solution, or is solvable by linear programming. Table 1 summarizes this information.

Set of restrictions (11) (14)

x? for Problem (8) emerges from closed-form expression (12). solving linear program (15) where right-hand side scalarPv? satisfies   (1− α1ˆ ) n∈N rn xˆn Λ b, n ≤ v? ≤ , 1 (1−p) αˆ

vector xˆ solves q linear program (19) and αˆ = 1 +

1+

1 1−p .

Table 1: Easily-solvable special cases of Problem (8)

AC

CE

PT

ED

Special case label 1 2

It is worth remembering that although both special cases may be quite restrictive in the con-

text of our motivating application, they provide useful insights for the design of solution methods addressing problem instances outside their scope. In the next section, we advance an exact solution algorithm for special case 2.

20

ACCEPTED MANUSCRIPT

3. An exact solution algorithm for the second special case n o On the assumption that we meet the requirements of special case 2 (see Table 1), let x[q]

0≤q≤Q

be a sequence of output vectors for some algorithm solving Problem (15) (i.e., x? = x[Q] ) where −1 Q ≥ 0. Because FP

0≤q≤Q

−1 other and guarantee the validity of v[q] LB ≤ F P

CR IP T

  (p) is monotonically increasing and bounded over p ∈ (0, 1) it e [q] ( ) θn xn n o n o follows that there must exist two bounded sequences v[q] , v[q] that approach each LB UB n∈N>0 x[q]

0≤q≤Q

(p) ≤ v[q] U B for all 0 ≤ q ≤ Q. However, e ? n∈N>0 (x? ) θn ( xn ) h i −1 (p), the algorithm can select some value v[q] ∈ v[q] in order to determine FP , v[q] ? LB U B e n∈N>0 (x? ) θn ( xn ) n o for each q ∈ {0, . . . , Q} such that the sequence v[q] encodes a family of approximations for 0≤q≤Q

Problem (15) each one being

AN US

X   z x[q] , min cn xn x

(31)

n∈N

subject to

(32)

xn ≥ ln for all n ∈ N.

(33)

M

h (x) ≥ v[q]

ED

To ascertain the feasibility and optimality status of x[q] resulting from Problem (31), the     [q]  h x[q] , and then compare the result to algorithm would use (30) to compute FP e θ x n n n∈N>0 (x[q] )    [q]  [q]  h x[q] p. If both FP ≥ p and v < v[q] e LB U B are valid, we would be guaranteed the θ x n∈N>0 (x[q] ) n n feasibility (albeit sub-optimality) of x[q] for Problem (8). Consequently, the next upper bound

CE

PT

[q+1] [q] v[q+1] U B would have to satisfy vU B < vU B in order to move the algorithm forward.     [q]  h x[q] < p would mean that x[q] is infeaBy contrast, the observation that FP e θ x n∈N>0 (x[q] ) n n sible for Problem (8), thereby requiring that v[q+1] > v[q] LB LB . Thanks to the monotone convergence n [q] o n [q] o [q] theorem for the pair of sequences vLB , vU B , the inequalities −∞ < v[q] ≤ LB ≤ v 0≤q≤Q

0≤q≤Q

v[q] U B < ∞ for q = 0, . . . , Q, and Problem (15)’s guarantee of optimality (when p < 1 −

rn kn ),

it is

AC

the case that the terminal iterative index Q must satisfy the conditions             

Q<∞ [q] −1 [Q] v[Q] = v[Q] = v? = F P LB = v U B = lim v q→Q

(34) (p) > 0. e θ (x? )

n∈N>0 (x? ) n

n

In other words, Problem (15) (which yields the optimal solution x? for Problem (8) in 21

ACCEPTED MANUSCRIPT

special case 2), effectively is the limiting version of Problem (31) as v[q] approaches v? = [Q] (p) where q = 0, . . . , Q and v? = v[Q] = v[Q] LB = v U B . Propositions 3 and 4 jointly   [0]   1 P (1− αˆ ) n∈N rn xˆn meaning that we can initialize v[0] state that Λ b, n ≤ v? ≤ 1 LB , Λ b, n , vU B , −1 FP

e ? n∈N>0 (x? ) θn (xn )

=

(1− α1ˆ )

(1−p) αˆ

P

n∈N rn xˆn 1 (1−p) αˆ

. The preceding insights underlie Algorithm 1 below.

CR IP T

v?U B

Algorithm 1 Exact solution algorithm for Problem (8) in special case 2 i.e., when p < 1 − 1: 2:

Initialize q := 0 and select some Φ ∈ (0, 1) [q] Initialize vcurr LB := vLB := Λ b, n

Solve Problem (19) to obtain xˆ , then set αˆ := 1 + P (1− α1ˆ ) n∈N rn xˆn [q] 4: Initialize vcurr := v := 1 UB UB

9: 10: 11: 12: 13: 14: 15: 16: 17: 18: 19: 20:

Initialize v := repeat Solve Problem (31) using v [q] to obtain x[q] , then store xcurr := x[q]  [q]   P [q] Evaluate F h x using (6), (30) e θ x n∈N>0 (x[q] ) n n     [q]  h x[q] ≥ p then if FP e θ x n∈N>0 (x[q] ) n n //x[q] is feasible for Problem (8) [q] [q] [q] Store vcurr LB := vLB , and update vU B := v [q] [q] [q] Decrease v by setting v := ΦvLB + (1 − Φ) v[q] else //x[q] is infeasible for Problem (8) [q] [q] [q] Store vcurr U B := vU B , and update vLB := v [q] [q] [q] Increase v by setting v := Φv + (1 − Φ) v[q] UB end if Increment q := q + 1 curr until vcurr LB = vU B
P Optimal solution x? := xcurr with cost z x? := n∈N cn xncurr

AN US

8:

(1−p) αˆ

+ (1 − Φ) v[q] UB

M

7:

Φv[q] LB

1 1−p

ED

5: 6:

[q]

q 1+

PT

3:

rn kn

CE

Lines 1 - 5 of Algorithm 1 use the results of Propositions 3 and 4 to initialize the algorithm’s

parameters, whereas lines 6 - 19 constitute a loop that at each iteration, determines the feasibility

AC

status of x[q] so as to either increase, or decrease v[q] (and simultaneously update its bounds) in preparation for the next oracle call to Problem (31). The validity of the condition in line 19 heralds the algorithm’s completion after a finite number of steps at which point the bounds on v[q] will have met each other (i.e., (34) will be fulfilled). Line 20 returns the optimal solution,

but if Φ ∈ (0, 1) in line 1 is chosen poorly, the convergence might be slower than is necessary, and so in practice, it would be sufficient to terminate the computation as soon as q satisfies some 22

ACCEPTED MANUSCRIPT

criterion such as

[q]

[q]

vU B −vLB

≤  stop for some small 0 <  stop  1.

[q] vU B

4. Complementary solution approaches

CR IP T

4.1. A gradient-based heuristic Our solution techniques thus far only guarantee optimality when p < 1 −

are only practically useful for the firm as long as 1 −

rn kn

rn kn .

Hence, they

is above, e.g., 0.95. When 1 −

rn kn

is not

high enough, the next best course of action is to seek low-cost feasible solutions for situations in r

which p ≥ 1 − knn . We propose Algorithm 2 as a heuristic procedure whose goal is to achieve that objective.

AN US

Algorithm 2 Heuristic procedure for Problem (8) in cases beyond Algorithm 1’s reach i.e., when r p ≥ 1 − knn r

r

Initialize π := 1 − knn −  where 0 <   1 − knn 2: Use Algorithm 1 to obtain an initial infeasible  solution x(π) :=      X  cn xn : FPn∈N (x) eθn (xn ) (h (x)) ≥ π, xn ≥ ln for all n ∈ N  arg min  with components    x  >0 n∈N x1 (π), . . . , x|N| (π)   cm cm c|N| =⇒ m1 = n 3: Construct the ordered set S , m1 , m2 , . . . , m|N| : r 1 ≤ r 2 ≤ . . . ≤ r m1 m2 |N| 4: Initialize i := 2, π previous := 0 and select some incentive step size ∆x > 0 (h (x(π))) 5: Use (6), (30) to evaluate F P e n∈N>0 (x(π)) θn (xn (π)) (h (x(π))) 6: while F P < p do e n∈N>0 (x(π)) θn (xn (π)) 7: Update π := FPn∈N (x(π)) eθn (xn (π)) (h (x(π)))

ED

M

1:

>0

−

c ∆x

mi cn log(1−b) ! rn 1 (1−p)3 rn log2 · kn 1−p

≥0

if π − π previous ≥ ∆π then Update xmi (π) := xmi (π) + ∆x else if i ≤ |N| − 1 then i := i + 1 end if Update xmi (π) := xmi (π) + ∆x end if Update π previous := π Use (6), (30) to evaluate FPn∈N (x(π)) eθn (xn (π)) (h (x(π))) >0 end while   P Feasible heuristic solution x?FHS := x(π) with cost z x?FHS := n∈N cn xn (π)

CE

9: 10: 11: 12: 13: 14:

Compute ∆π :=

PT

8:

AC

15: 16: 17: 18:

19: 20:

23

ACCEPTED MANUSCRIPT

Algorithm 2 works as follows. Choose an appropriate initial minimum reliability π and use it to define a relaxed problem that falls under the jurisdiction of special case 2 (lines 1 - 2). The optimal solution x(π) for that relaxed problem will certainly be infeasible for Problem (8) thanks to π < p where p is the minimum reliability that we actually require, and p ≥ 1 −

rn kn .

Then, in a

CR IP T

new set S, sort the classes of N such that they are in ascending order of their cost-responsivity

ratios (line 3). Next, initialize some index i := 2, some probability value π previous := 0, choose

some strictly positive step size ∆x, and then use (6) and (30) to evaluate the chance constraint’s left-hand side FPn∈N

e

θ (x (π)) >0 (x(π)) n n

As long as FPn∈N

(h (x(π))) at the incumbent solution x(π) (lines 4 - 5).

e

θ (x (π)) >0 (x(π)) n n

(h (x(π))) is strictly below p, remain in the loop beginning on

line 6 and ending on line 19. Once inside that loop, update π by assigning it the value of e

θ (x (π)) >0 (x(π)) n n

(h (x(π))) then use the incumbent index i to compute some non-negative value

AN US

FPn∈N

∆π (lines 7 - 8). Referring to (13) and line 8 of Algorithm 2, we observe that ∆π’s denominator is the chance constraint’s shadow price on the assumption that (11) holds. Additionally, ∆π’s numerator is the extra cost for increasing the incentive afforded to class mi ∈ S by ∆x units. Since (11) does not hold under the circumstances (p ≥ 1 −

rn kn )

for which Algorithm 2 is designed, it

M

stands to reason that ∆π represents only an estimate of the fair (i.e., expected) increase in the chance constraint’s left-hand side when we afford ∆x extra incentives to class mi . Contrastingly, π − π previous is the observed increase in the chance constraint’s left-hand side due to the provision

ED

of ∆x new incentives to class mi .

As a consequence, Algorithm 2 interprets the validity of π − π previous ≥ ∆π as sufficient

PT

evidence for the hypothesis that class mi is attractive enough to continue incentivizing. If π −

π previous ≥ ∆π, then update the incumbent solution x(π) by adding ∆x to class mi ’s incentive (lines 9 - 10). However, if π − π previous < ∆π, verify that the incumbent class mi is not the class

CE

m|N| with the greatest cost-responsivity ratio. If that verification yields a positive result, move on

to the next class in S (lines 11 - 14). Furthermore, (still under the “if π − π previous < ∆π” branch)

AC

update the incumbent solution x(π) by adding ∆x to class mi ’s incentive (line 15).

Finally, update π previous by assigning it π’s current value, then use (6) as well as (30) to

evaluate the chance constraint’s left-hand side FPn∈N

solution x(π) (lines 17 - 18). If FPn∈N

e

θ (x (π)) >0 (x(π)) n n

e

θ (x (π)) >0 (x(π)) n n

(h (x(π))) at the new incumbent

(h (x(π))) ≥ p, exit the loop because the new

incumbent solution x(π) is no longer infeasible for Problem (8) (line 20). In short, Algorithm 2’s operating philosophy is to start from an infeasible point, and then navigate along the cheapest 24

ACCEPTED MANUSCRIPT

computable route to the feasible space by iteratively updating a function of the chance constraint’s estimated shadow price gradient. Another approach to our CCNLP when 1 −

rn kn

is not nearly as high as we want it to be, is

to adopt the sample average approximation perspective of Luedtke and Ahmed (2008). In that

CR IP T

asymptotically optimal approach, we replace the chance constraint with a (ideally) large set of

deterministic constraints featuring empirical realizations of the CCNLP’s random parameters. The next subsection explains this method in detail.

X   cn xn z x? = min x

n∈N

AN US

4.2. An asymptotically optimal mixed integer linear program ! Y X 1 (e , we may reformulate By virtue of 1 − an )−rn xn ≥ b ⇐⇒ rn xn log (e an ) ≥ log 1−b n∈N n∈N  ! X 1    ≥ p and consequently Problem (1) to (3) as an ) ≥ log rn xn log (e (2) as Pr  1−b  n∈N

M

subject to  ! X  1  ≥ p an ) ≥ log Pr  rn xn log (e 1−b  n∈N

ED

xn ≥ ln for all n ∈ N.

(35)

(36) (37)

PT

  1 Note that b ∈ (0, 1) =⇒ log 1−b > 0 and recall that e an ∈ A is uniformly-distributed on  the continuous support 0, a+n where a+n > 1. Now for every n ∈ N, suppose that a[d] n is the dth

empirical realization of e an in a Monte Carlo sample of size D ≥ 1. Using the sample average

AC

CE

approximation approach of Luedtke and Ahmed (2008), we obtain for (36)’s left-hand-side    ! Y   X 1  −rn xn     (e Pr 1 − an ) ≥ b = Pr  rn xn log (e an ) ≥ log 1−b  n∈N n∈N  ! D   1 X X 1  [d]  = lim I  rn xn log an ≥ log D↑∞ D 1−b  d=1 n∈N

(38)

where I(·) is the indicator function taking on the value of one when its argument is true, and ! X 1    zero otherwise. Alternatively stated, (36)’s left-hand-side is the average of I  rn xn log (e an ) ≥ log 1−b  n∈N
taken over infinitely many empirical realizations of the set A = {e an }n∈N . Consequently, z x? = 25

ACCEPTED MANUSCRIPT

    lim e zS A D, x?S A such that e zS A D, x?S A is the value of the following mixed integer linear program

D↑∞

(MILP) with D binary variables, |N| continuous variables, and 1 + D + |N| constraints. X   e zS A D, x?S A , τmin cn xn ,...,τ x

D

(39)

n∈N

subject to D

1 X τd ≥ p D d=1

CR IP T

1

(40)

    X   I ∃n ∈ N : log a[d] >0 · rn xn log a[d] ≥ τd log n n

(41)

AN US

n∈N

! 1 for d = 1, . . . , D 1−b

τd ∈ {0, 1} for d = 1, . . . , D

(42)

xn ≥ ln for all n ∈ N.

(43)

The expression in (38) gives rise to (40) and (41), and as the sample size D approaches infinity, Problem (39) better approximates Problem (1). We need the indicator multiplier on the

M

left-hand side of (41) to avoid being compelled to set all x decision variables to zero when there n  o exists a d ∈ {1, . . . , D} for which all elements of the realized set log a[d] are non-positive. n n∈N

ED

Addressing this technicality is essential because then if at least one element of {ln }n∈N is strictly

positive, the imperative to force all x variables to zero will lead to Problem (39) being infeasible.

PT

5. Numerical experiments and discussion Figure 1 is a concise illustration of the ensemble of methods that we propose for tackling

AC

CE

Problem (1).

26

ACCEPTED MANUSCRIPT

Problem (1) is a CCNLP with optimal solution x? , optimal cost z x? and p ∈ (0, 1) as the right-hand side of its chance constraint.




Problem (1) is equivalent to Problem (8) thanks to (7). )

CR IP T

:=

cn arg min n∈N rn

Algorithm 2 yields x?FHS which is always feasible for Problem (8). Alternatively, for each n ∈ N, generate a Monte n Carlo o No (Branch B) sample a[d] n 1≤d≤D , then solve Problem (39) to determine  e zS A D, x?S A where   zS A D, x?S A = lim e D↑∞   z x? . However, the realized x?S A is only potentially feasible for Problem (8).

AN US

Determine n

(

PT

ED

M

Algorithm 1 yields the optimal solution x? . r If ln = 0 for all n , n, p < 1 − knn ? then x? already exists Yes (Branch A) in the closed-form expressed in (12).

Figure 1: Options for handling Problem (1)

CE

We investigated a real problem instance for a firm with |N| = 9 customer classes that required

an aggregate product desirability index of at least b = 0.95. The firm assigned to each class n ∈ N

AC

an equal weight of kn =

1 |N|

= 19 , with cn , a+n , and ln listed for all n ∈ N in Table 2 below.

27

ACCEPTED MANUSCRIPT

Value ($ day) 550.1615 890.1414 620.3598 890.4996 1250.0736 330.5905 1080.5883 700.2470 1120.3008

Parameter a+1 a+2 a+3 a+4 a+5 a+6 a+7 a+8 a+9

Value 6.6094 8.3401 4.6046 7.3375 6.7474 3.6565 4.9438 2.2818 4.4864

Table 2: Real parameter settings

Parameter l1 l2 l3 l4 l5 l6 l7 l8 l9

Value (days) 0.4000 0.5500 0.4000 0.2000 0.3500 0.2000 0.5000 0.4500 0.3000

CR IP T

Parameter c1 c2 c3 c4 c5 c6 c7 c8 c9

AN US

We coded all algorithms in C++ and ran all experiments on a 3.4 GHz, 32 GB RAM computer. Also, we employed the CPLEX 12.6 commercial solver to solve any arising linear, or mixed integer linear programs.

5.1. Real problem instance: Branch A of Figure 1

r6 k6

= 1−

1 a+6

M

and 1 −

(

) ( + ) a cn cn = arg min n =6 n∈N rn n∈N kn = 0.7265. Hence, as long as p does not attain 0.7265, we can follow

The parameter settings in Table 2 give us the index n = arg min

Branch A in Figure 1 to determine the firm’s optimal solution. We implemented Algorithm 1

ED

using Φ := 0.618 which is the inverse golden ratio specified to three decimal places. Table 3 displays the cost and runtime performance for p = 0.05 to p = 0.70 < 0.7265 in incremental steps of 0.05.

PT


Optimal cost, z x? ($) 28,641.1 29,816.1 31,179.5 32,784.3 34,704.7 37,049.2 39,983.9

AC

CE

p 0.05 0.10 0.15 0.20 0.25 0.30 0.35

Runtime (secs) 0.187 0.141 0.140 0.125 0.124 0.156 0.109

p 0.40 0.45 0.50 0.55 0.60 0.65 0.70


Optimal cost, z x? ($) 43,775.6 48,882.8 56,167.5 67,465.4 87,507.7 133,378.0 350,956.0

Runtime (secs) 0.109 0.124 0.125 0.109 0.094 0.062 0.078

Table 3: Optimal results (Algorithm 1) for the real problem

Algorithm 1 computed x? extremely quickly which is unsurprising given that its main job is

to sequentially solve linear programs. The results in Table 3 confirm the intuition that raising the minimum desirable reliability p makes the chance constraint harder to satisfy, thereby increasing 28

ACCEPTED MANUSCRIPT

the optimal cost. That is, the price that the firm pays for limiting its infeasibility risk (i.e., being more conservative) is inexorably a higher minimum cost. Interestingly, the optimal cost rises fairly slowly from p = 0.05 up until about p = 0.60 at which point it shoots up rather steeply. In to the tradeoff between the optimal cost and p. 5.2. Real problem instance: Branch B of Figure 1

CR IP T

practice, what this means for the firm is that as long as p is low, it need not pay much attention

Since our exact solution techniques fail for values of p ≥ 0.7265 = 1 −

rn kn

= 1−

r6 k6 ,

we

ran Algorithm 2 for p = 0.75 to p = 0.95 in steps of 0.05 all along using an incentive step size

∆x = 5 and an initial offset  = 0.01. Table 4 displays the solution costs as well as the runtime

  Cost, z x?FHS ($) 1,002,730.2 1,887,700.8 2,387,780.3 4,144,540.1 13,337,800.0

Chance constraint’s left-hand side value at x?FHS 0.750058 0.800178 0.850004 0.900128 0.951501

Runtime (secs) 1.094 2.672 3.281 4.391 13.969

M

p 0.75 0.80 0.85 0.90 0.95

AN US

performance of Algorithm 2.

ED

Table 4: Heuristic results (Algorithm 2) for the real problem

To determine the values in Table 4’s third column, we invoked (6) and (30). As is evident from Table 4’s results, all instances for which we deployed Algorithm 2 produced a feasible solution

PT

x?FHS in under 14 seconds. To boot, the chance constraint’s surplus value at x?FHS (measured by the probability in Table 4’s third column minus that in its first column) in no case exceeded

CE

0.001501. This observation makes us reasonably confident that for p = 0.75 to p = 0.95,  
Algorithm 2’s cost z x?FHS might be close to that of that the (unknown) optimum z x? . As a counterweight to Algorithm 2, we solved the sample average approximation MILP

AC

(Problem (39)) using Monte Carlo samples of size D = 15, 000 for each random parameter e an . Table 5 summarizes Problem (39)’s results.

29

ACCEPTED MANUSCRIPT

  Cost, e zS A D, x?S A ($) 83,908.3 89,403.9 94,899.7 100,434.0 106,011.0

Chance constraint’s left-hand side value at x?S A 0.566776 0.583830 0.606811 0.616658 0.623950

Runtime (secs) 12,609.400 11,315.100 11,498.200 12,295.900 13,250.100

CR IP T

p 0.75 0.80 0.85 0.90 0.95

Table 5: sample average approximation (Problem (39)) results with D = 15, 000 for the real problem. The choice of D = 15, 000 was driven by pragmatic computational considerations. We began with D = 1000 and then incremented D in steps of 1000 until every instance of p ∈ {0.75, 0.80, 0.85, 0.90, 0.95} took more than 10, 800 seconds (i.e., three hours) to yield a solution (by D = 15, 000, each run essentially consumed all available RAM memory).

Unfortunately, solving the sample average approximation MILP took on average 12, 193.74

AN US

seconds and yet in each instance, the result was an infeasible solution that is quite distant from the

feasible region (note that the probabilities in the third column of Table 5 are always significantly lower than those in its first column). Surmisedly, the sample average approximation’s subpar performance stems from the fact that the entries in (41)’s left-hand side coefficient matrix

ED

M

      r1 log a[1] r2 log a[1] ... 1 2   [2]   [2]    r1 log a1 r2 log a2 ... Ξ|N| (D) ,  .. .. ..  . . .    [D]   [D]  r1 log a1 r2 log a2 ...

   r|N| log a[1]  |N|   [2]   r|N| log a|N|   ∈ RD×|N| ..   .  [D]  r|N| log a|N|

(44)

are scaled logarithms of random variable realizations (as opposed to mere linear transfor-

PT

mations of random variable realizations). This detail is relevant because the slowly-changing (but nonlinear) nature of the logarithm function means that the collective disparity within the ensemble of Ξ|N| (D)’s rows is likely to be quite small even for large values of the sample size

AC

CE

D. In other words, the logarithmic transformation could possibly have the implication that  ! D X X     1 1   toward in the limit as D −→ ∞, the convergence of I  rn xn log a[d] ≥ log n D d=1 n∈N 1−b     ! Y   X 1  −rn xn  (see (38), (41)) is orders of mag(e Pr 1 − an ) ≥ b = Pr  rn xn log (e an ) ≥ log 1−b  n∈N n∈N     D  X  1 X X [d]   I  nitude slower than the situation in e.g., lim rn xn an ≥ · = Pr  rn xne an ≥ ·. The D↑∞ D d=1 n∈N n∈N message therefore seems to be that for our application, D = 15, 000 is a tiny sample size in spite of the significant computational resources that it requires. In other words, our sample average approximation (Problem (39)) requires an enormous sample size to serve as a reasonable sur30

ACCEPTED MANUSCRIPT

rogate approximation for Problem (1). Overall, the experimentation for situations outside the special cases (i.e., for p = 0.75 to p = 0.95) strongly suggests that despite the sample average approximation approach’s asymptotic optimality, heuristic Algorithm 2 is the more pragmatic

CR IP T

option for managing the firm’s new product risk. 5.3. Artificial problem instances

To gauge the impact of each class’ maximum propensity a+n on the optimal cost, we retained     the c = c1 c2 . . . c9 , l = l1 l2 . . . l9 parameter vectors in Table 2 but then for b = 0.15 to b = 0.95 in steps of 0.20, and for p = 0.05 to p = 0.95 in steps of 0.05, we randomly generated new problem instances segregated into two regimes namely, “low propensity”, and

AN US

“high propensity”. Moreover, to create a domain of p that is large enough to accommodate

p = 0.95, we enforced an adjusting mechanism in the synthetic data generation process to ensure settings.

rn kn

≤ 0.99 was valid for each problem instance. Table 6 below shows the various regime

Propensity regime Low propensity (LP) High propensity (HP)

 Each member of a+n n∈N is independently drawn from uniform distribution supported on (1, 5] uniform distribution supported on (10, 15]

M

that 1 −

ED

Table 6: Maximum propensity regimes for artificial problem instances

In the low propensity regime, for each combination of b and p, we generated 200 independent

PT

problem instance replications. In each replication, we independently drew every member of the  set a+n n∈N from a uniform distribution supported on (1, 5]. The performance measure for each

CE

instance in the low propensity regime was the instance’s sample average cost which we denote 1 P200 ?
by 200 rep=1 zrep,LP x . For the high propensity regime, we followed the same approach as

AC

before but instead drew each a+n independently, and uniformly from (10, 15]. Then, we recorded 1 P200 ?
the sample average cost 200 of each high propensity instance. We summarize rep=1 zrep,HP x the results using the following table.

31

ACCEPTED MANUSCRIPT

b = 0.35 18.65234275 28.54536050 31.45232903 39.84527022 47.41604323 48.34189085 49.89741580 50.04284784 52.22787066 53.20446036 56.84161336 61.98109535 63.84824077 59.72515046 63.46285344 65.22197593 68.73934021 83.47208062 99.56628709

1 200 1 200

P200

rep=1 zrep,LP

P200

(x? )

rep=1 zrep,HP

b = 0.55 22.78819931 23.98395825 27.76585338 35.47892728 31.70049550 41.27825548 33.16643303 31.73222866 40.32051815 37.23292449 36.44048319 40.37320402 48.71982145 46.34528347 43.62078402 44.05877410 42.78931223 53.10911033 60.78208297

(x? )



b = 0.75 18.33424511 21.18484231 20.07511457 23.28882448 22.84617013 25.66746361 21.66468357 29.62269987 31.84347743 26.13559732 28.06066021 27.95806517 27.45823559 29.31122299 29.21574962 29.94583538 33.70967891 32.90040030 36.43561013

AN US

0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95

b = 0.15 8.17901248 16.71791714 31.46531821 36.89015512 27.55611659 50.42778259 49.46485215 65.20290175 61.76805981 70.39703950 87.64302878 83.93771269 81.04582869 100.05862195 102.76844043 109.02134009 126.87541709 133.11756672 209.52428472



b = 0.95 9.63796066 12.97998315 16.45993035 13.58866495 13.93057354 14.58038868 14.37653667 15.38164002 13.45553917 13.52168486 14.06734829 15.59069355 15.80561502 12.92743419 16.02836000 15.31770280 16.78372742 18.16681308 16.49887672

CR IP T

1000 log

p

M

Table 7: 1000log-ratio of sample average optimal costs for combinations of b and p (artificial problem instances)

200

ED

In a nutshell, Table 7 suggests that b strongly controls the impact of p on the ratio  1 P200  zrep,LP (x? ) log 2001 Prep=1 which is a comparison of the average costs of the different propensity 200 (x? ) z rep=1 rep,HP

regimes. When b is quite high, the comparative ratio of average optimal costs is relatively insen-

PT

sitive to p. In other words, when the firm requires a high minimum aggregate product desirability index, the reliability parameter is relatively inconsequential to the impact of customers’ propensity on the firm’s optimal cost. As the minimum aggregate product desirability index decreases,

CE

the optimal cost impact of customers’ propensity tends to react more sensitively to the reliability parameter. Additionally, all values in Table 7 are positive. That is,

AC

  1000 log 

1 200 1 200

P200

rep=1 zrep,LP P200 rep=1 zrep,HP


 200 200 X     x?  1 X  > 0 ⇐⇒ 1 zrep,LP x? > zrep,HP x? 200 rep=1 200 rep=1 (x? )

for all values of b and p in Table 7. Unsurprisingly, it appears therefore that the firm’s optimal

incentivization cost is inversely related to its customers’ propensity.

32

ACCEPTED MANUSCRIPT

6. Conclusion In this paper, we adopted a chance-constrained optimization framework to address a nonlinear program plagued by multiplicative noise. Chance-constrained models are a subset of non-

CR IP T

convex programs which in and of themselves constitute a class of quite difficult non-traditional optimization problems. We proved the existence of a weakly restrictive special case for which it is possible to determine the optimal solution by solving a sequence of appropriately-designed

linear programs. In this regard, we showed that tackling a certain category of nonlinear programs with random parameters is conditionally no more difficult than solving several (easy) linear op-

timization problems. We then designed a solution algorithm for the aforementioned special case

AN US

and applied it to the real problem of managing a consumer firm’s new product marketing risk. In a strongly restrictive special case, we were able to express our problem’s solution in closed-form.

For problem instances outside the arena of the special cases, we used insights from the theoretical analysis to advance a heuristic which we compared to a Monte-Carlo-based sample average approximation method (known to be optimal in the limit of an infinite sample size). Computational experiments using the firm’s data yielded strong evidence for the practical superiority of

M

our heuristic over the sample average approximation approach. To conclude the investigation, we created and solved synthetic problem instances (covered by the weakly restrictive special case)

ED

so as to observe the impact of varying a few model parameters. We believe that there are at least three worthwhile extensions of this research. First, broadening the analysis to include scenarios of non-zero mutual correlation between the random elements

PT

in the chance constraint would serve to greatly increase our model’s scope of usefulness. Second, the static nature of our chance-constrained program means that there is no interactive feedback

CE

between the firm and its customers. A significant future contribution would be to model the firm-customer relationship in the framework of closed-loop dynamic programming, and to determine cost-effective policies for providing product promotion incentives based on observed state

AC

information. Finally, the fact that we can−for some special cases−solve our problem via linear programming (which is a special case of convex programming), suggests that in general, perhaps it might be possible to solve our problem via a convex programming algorithm. We leave the resolution of this last conjecture to future research.

33

ACCEPTED MANUSCRIPT

Acknowledgments The authors are grateful to both anonymous referees for their constructive and insightful comments which significantly improved the paper’s quality. This research was partially supported by

CR IP T

the University of Mary Washington Faculty Development Fund. References References

Ahmed S., Shapiro A. 2008. Solving chance-constrained stochastic programs via sampling and integer programming. State-of-the-Art Decision-Making Tools in the Information-Intensive Age, Chapter 12, pp. 261-269

AN US

Araya A., Kisekka I., Gowda P. H., Prasad P. V. V., 2018. Soybean crop-water production functions in a humid region across years and soils determined with APEX model. Agricultural Water Management 204 180–191 Bairam E., 1994. Homogeneous and nonhomogeneous production functions: theory and applications Avebury; Brookfield, USA: Ashgate, Aldershot.

Blackmore L., Ono M., Williams B. C., 2011. Chance-constrained optimal path planning with obstacles. IEEE Transactions on Robotics 27(6) 1080–1094

Charnes A., Cooper W. W., Symmonds, G. H., 1958. Cost horizons and certainty equivalence: an approach to stochastic

M

programming of heating oil. Management Science 4 235–263

Charnes A., Cooper W. W., 1959. Chance-constrained programming. Management Science 6(1) 73–79 Dentcheva D., Pr´ekopa A., Ruszczynski A., 2000. Concavity and efficient points of discrete distributions in probabilistic

ED

programming. Mathematical Programming 89(1) 55–77 Doppa J. R., Yu J., Tadepalli P., Getoor L., 2010. Learning algorithms for link prediction based on chance constraints. Proceedings of the 2010 European Conference on Machine Learning and Knowledge Discovery in Databases: Part I,

PT

Springer-Verlag, Berlin, Heidelberg

Fakhar M., Mahyarinia M. R., Zafarani J., 2018. On nonsmooth robust multiobjective optimization under generalized convexity with applications to portfolio optimization. European Journal of Operational Research 265(1) 39–48

CE

Humphrey T. M., 1997. Algebraic production functions and their uses before Cobb-Douglas. FRB Richmond Economic Quarterly 83(1) 51–83 Available at SSRN: https://ssrn.com/abstract=2129863

Geletu A., Kl¨oppel M., Zhang H., Li P., 2013. Advances and applications of chance-constrained approaches to systems optimisation under uncertainty. International Journal of Systems Science 44(7) 1209–1232

AC

Grassetti F., Hunanyan G., 2018. On the economic growth theory with Kadiyala production function. Communications in Nonlinear Science and Numerical Simulation 58 220–232

Kargarian A., Fu Y., Wu H., 2016. Chance-constrained system of systems based operation of power systems. IEEE Transactions on Power Systems 31(5) 3404–3413

Li P., Arellano-Garcia H., Wozny G., 2008. Chance constrained programming approach to process optimization under uncertainty. Computers & Chemical Engineering 32(1) 25–45

34

ACCEPTED MANUSCRIPT

Luedtke J., Ahmed S., 2008. A sample approximation approach for optimization with probabilistic constraints. SIAM Journal on Optimization 19(2) 674–699 Luedtke J., Ahmed S., Nemhauser G. L., 2010. An integer programming approach for linear programs with probabilistic constraints. Mathematical Programming 122(2) 247–272 Available at SSRN: https://ssrn.com/abstract=1749083

CR IP T

Mishra S. K., 2010. A brief history of production functions. The IUP Journal of Managerial Economics 3(4) 6–34 Nemirovski A., Shapiro A., 2006. Scenario Approximations of Chance Constraints. Probabilistic and Randomized Methods for Design under Uncertainty, Springer, London

Nemirovski A., Shapiro A., 2007. Convex approximations of chance constrained programs. SIAM Journal on Optimization 17(4) 969–996

Pagnoncelli B. K., Ahmed S., Shapiro A., 2009. Sample average approximation method for chance constrained programming: theory and applications. Journal of Optimization Theory and Applications 142 399–416

AN US

Pr´ekopa A., 1973. Contributions to the theory of stochastic programming. Mathematical Programming 4(1) 202–221

Pr´ekopa A., 1993. Programming under probabilistic constraint and maximizing a probability under constraints. Rutcor Research Report, RRR 35-93, Rutgers University, Brunswick, NJ, USA

Pr´ekopa A., 1995. Stochastic Programming Dordrecht, The Netherlands: Kluwer Academic Publishers. Pr´ekopa A., Vizv´ari B., Badics T., 1998. Programming under probabilistic constraint with discrete random variable. New Trends in Mathematical Programming, eds. S. Vajda, F. Giannessi, and T. Rapcsak, Dordrecht, The Netherlands: Kluwer, pp. 235–255 Operations Research Letters 39 188–192

M

Pr´ekopa A., Yoda K., Subasi M. M., 2011. Uniform quasi-concavity in probabilistic constrained stochastic programming. Sepp¨al¨a Y., 1972. A chance-constrained programming algorithm. BIT Numerical Mathematics 12(3) 376–399

ED

Shahbaz M., Benkraiem R., Miloudi A., Lahiani A., 2017. Production function with electricity consumption and policy implications in Portugal. Energy Policy 110 588–599 Wei Y., Bevers M., Belval E., Bird B., 2015. A chance-constrained programming model to allocate wildfire initial attack resources for a fire season. Forest Science 61(2) 278–288

PT

Wendt M., Li P., Wozny G., 2002. Nonlinear chance-constrained process optimization under uncertainty. Industrial & Engineering Chemistry Research 41(15) 3621–3629 Wibe S., 1984. Engineering production functions: a survey. Economica 51(204) 401–411

CE

Zhang H., Li P., 2011. Chance constrained programming for optimal power flow under uncertainty. IEEE Transactions on Power Systems 26(4) 2417–2424

Zhang Y., Jiang R., Shen S., 2016. Ambiguous chance-constrained bin packing under mean-covariance information.

AC

arXiv:1610.00035v2

Zhang B., Feng G., Ahuja L. R., Kong X., Ouyang Y., Adeli A., Jenkins J. N., 2018. Soybean crop-water production functions in a humid region across years and soils determined with APEX model. Agricultural Water Management 204 180–191

35

ACCEPTED MANUSCRIPT

Appendix A. Proof of Proposition 1

X

! e an −rn xn log + . an (x)

X

n∈N>0

AC

CE

PT

ED

M

AN US

n∈N>0 (x)

e θn (xn ) =

CR IP T

If e θn (xn ) ∈ T (x) is exponentially-distributed with rate rn1xn , then by the probability integral    
e e transform, the random variable e an , 0 + a+n − 0 exp − θrnn(xxnn ) = a+n exp − θrnn(xxnn ) is uniformly    e distributed over the domain 0, a+n . Rearranging a+n exp − θrnn(xxnn ) to make e θn (xn ) the subject of the  ea  equation gives us e θn (xn ) = −rn xn log a+nn and because of the mutual independence assumption h i for all elements of A = {e an }n∈N , it follows that Cov e θn (xn ) , e θm (xm ) = 0 for all n , m and

36

(A1)

ACCEPTED MANUSCRIPT

Consequently,

n∈N

AC

CE

PT

ED

M

AN US

CR IP T

    X e θn (xn ) ≤ h (x) ≥ p FPn∈N (x) eθn (xn ) (h (x)) ≥ p ⇐⇒ Pr  >0 n∈N>0 (x)  ! X  X  r n e  ≥ p by (6) θn (xn ) ≤ log(1 − b) − rn xn log ⇐⇒ Pr  kn  n∈N>0 (x) n∈N>0 (x)  ! X  X  1 kn  e θn (xn ) ≤ log(1 − b) − rn xn log +  ≥ p since rn = + ⇐⇒ Pr  an an n∈N>0 (x) n∈N>0 (x)   X  X
 e θn (xn ) ≤ log (1 − b) + rn xn log a+n  ≥ p ⇐⇒ Pr  n∈N>0 (x) n∈N>0 (x)   X X  X

 + +  e θn (xn ) ≤ log (1 − b) − log an + (rn xn + 1) log an  ≥ p ⇐⇒ Pr  n∈N>0 (x) n∈N>0 (x) n∈N>0 (x)   ! X  X 
1 − b + e ⇐⇒ Pr  θn (xn ) ≤ log Q (rn xn + 1) log an  ≥ p + + n∈N>0 (x) an n∈N>0 (x) n∈N>0 (x)   ! ! X  X
−(rn xn +1)  e an 1−b +  ≥ p by (A1) log an − ⇐⇒ Pr  −rn xn log + ≤ log Q + an n∈N>0 (x) an n∈N>0 (x) n∈N>0 (x)  ! !−r x X  X 
−(rn xn +1) e an n n 1−b + ⇐⇒ Pr  log + log an ≤ log Q +  ≥ p + an n∈N>0 (x) an n∈N>0 (x) n∈N>0 (x)   !−r x   Y e
−(rn xn +1) an n n Y 1−b +  ≥ p  Q ⇐⇒ Pr  a ≤ · n  + + a a (x) n∈N n n >0 n∈N>0 (x) n∈N>0 (x)   Q ! Y   an )−rn xn 1−b n∈N>0 (x) (e −r x n n 1 −  ≥ p (e ) Q Q ⇐⇒ Pr ≤ a ≥ p ⇐⇒ Pr ≥ b n   + + n∈N>0 (x) an n∈N>0 (x) an n∈N>0 (x)   Y   (e ⇐⇒ Pr 1 − an )−rn xn ≥ b ≥ p since (e an )0 = 1 for all n ∈ N.

37

ACCEPTED MANUSCRIPT

Appendix B. Proof of Proposition 2 ( ) cn be such that Let n = arg min n∈N rn if n = n

(B1)

CR IP T

           log(1−b)     > 0 l max ,   n rn   1   rn log kn · 1−p   ? ∞ > xn =       0

otherwise.

Additionally, assume the veracity of (11). It suffices to show that the assumption that (11) is X   e true leaves us no choice but to conclude that (B1) is optimal for Problem (8) when θn xn? = n∈N>0 (x? )

1 rn xn?

and possesses the

AN US

    e θn xn? . Note that e θn xn? is exponentially-distributed with rate parameter following quantile function that is linear with respect to

xn?

Feθ−1 x? (p) = −rn xn? log(1 − p). n( n ) By (B1), we obtain

(B2)

PT

ED

M

           log (1 − b)    ?   ∞ > xn = max  l ,  n     rn 1      rn log kn · 1−p  ! ! rn rn ? ⇐⇒ log (1 − b) = rn xn log − log(1 − p) for xn? ≥ ln , and p < 1 − kn kn ! rn rn ? ⇐⇒ log (1 − b) − rn xn log = −rn xn? log(1 − p) for xn? ≥ ln , and p < 1 − kn kn   rn ⇐⇒ h xn? = Feθ−1 x? (p) for xn? ≥ ln , and p < 1 − by (B2) and (6). n( n ) kn

) cn But p < 1 − is one of the conditions in (11), and by virtue of n = arg min and (7), n∈N rn it follows that (B3) is both a necessary and sufficient condition for the following optimization

CE

(

(B3)

AC

rn kn

38

ACCEPTED MANUSCRIPT

problem: X   cn xn z x? = min x

(B4)

n∈N

CR IP T

subject to Feθn ( xn ) (h (x)) ≥ p ⇐⇒ h (x) ≥ Feθ−1 x (p) n( n)

(B5)

xn ≥ ln = 0 for all n , n, where n ∈ N.

(B7)

xn ≥ ln

(B6)

AN US

We observe that Problem (B4) (which is linear thanks to the linearity of h (x), and Fe−1 (p) = θn ( xn ) −rn xn log(1 − p) in the x decision variables) has its constraint sets (B5) and (B7) totally depen-

dent on the validity of (11). Besides, Problem (B4) is just a special case of Problem (8) given X   e that θn (xn ) = e θn xn . Hence, the assumption that (11) is true leaves us no choice but to n∈N>0 (x) X     e conclude that (B1) is optimal for Problems (8) and (B4) when θn xn? = e θn xn? .

M

n∈N>0 (x? )

Appendix C. Proof of Proposition 3

dominance gives us

ED

For any specific n ∈ N note that the relationship between variance and first-order stochastic

PT

   X  h  i  e Var  θn (xn ) ≥ Var e θn xn =⇒ FPn∈N (x) eθn (xn ) (h (x)) ≤ Feθn ( xn ) (h (x)) >0

CE

n∈N>0 (x)

−1 where FP

e

AC

n∈N>0 (x) θn (xn )

−1 =⇒ FP

e

n∈N>0 (x) θn (xn )

(p) ≥ Feθ−1 x (p) n( n)

(C1)

(0) = Fe−1 (0) = 0. Now, consider the following reformulation of θn ( xn )

39

ACCEPTED MANUSCRIPT

Problem (8): X   cn xn z x? = min x,v

(C2)

n∈N

−1 h (x) ≥ v = FP

e

n∈N>0 (x) θn (xn )

(p)

xn ≥ ln for all n ∈ N.

CR IP T

subject to (C3)

(C4)

Combining (C1), (C3), and (B2) yields the following equivalent form of Problem (C2)

x,v

AN US

X   z x? = min cn xn

(C5)

n∈N

subject to

rn kn

⇐⇒ log

feasible for (C6) because



rn kn

·

(C6)

xn ≥ ln for all n ∈ N.

(C7)

1 1−p



< 0, it turns out that (14) ensures that at least one x is

M

Since p < 1 −

h (x) ≥ v ≥ −rn xn log(1 − p)

AC

CE

PT

ED

! !   rn rn rn 1 1 p<1− =⇒ ∃xn ∈ [0, ∞) : Λ b, n ≥ rn xn log · thanks to log · <0 kn kn 1 − p kn 1 − p ! ! X rn rn − rn xn log ≥ −rn xn log(1 − p) =⇒ ∃xn ∈ [0, ∞) : log(1 − b) − rn ln log kn kn n∈N\n  ! !    X   rn rn    log(1 − b) − r x log − r x log ≥ −rn xn log(1 − p) =⇒ ∃xn ∈ [0, ∞) : inf   n n n n    xn ≥ln kn  kn   (x)\n n∈N >0 n∈N>0 (x)\n ! X rn =⇒ ∃xn ∈ [0, ∞) : inf log(1 − b) − ≥ −rn xn log(1 − p) rn xn log xn ≥ln kn n∈N (x) =⇒ ∃xn ∈ [0, ∞) :

n∈N>0 (x)\n

inf

xn ≥ln n∈N>0 (x)\n

>0

h (x) ≥ −rn xn log(1 − p)

=⇒ ∃xn ∈ [0, ∞) : h (x) ≥ −rn xn log(1 − p) thereby ensuring the validity of (C6). (

) cn guarantees the existence As a result, (14) in conjunction with the definition n = arg min n∈N rn 40

ACCEPTED MANUSCRIPT

of a v? that is optimal for Problem (C5) and satisfies (C6) such that  ! !     X  rn r   n  ? log(1 − b) − − r r x log x log ≥ v? ≥ −rn xn? log(1 − p) ∞ > inf h x? = inf   n n n n     xn ≥ln xn ≥ln k k   n n n∈N\n n∈N>0 (x)\n n∈N>0 (x)\n ! ! X rn rn ? ∞ > log(1 − b) − rn ln log − rn xn log ≥ v? ≥ −rn xn? log(1 − p) k k n n n∈N\n !   rn =⇒ ∞ > Λ b, n − rn xn? log ≥ v? ≥ −rn xn? log(1 − p) kn   =⇒ Λ b, n ≤ v? < ∞ (C8) 

CR IP T



AN US

Therefore, we may rewrite Problem (C5) as X   z x? = min cn xn x

(C9)

n∈N

subject to

h (x) ≥ v?

M

xn ≥ ln for all n ∈ N

But Problem (15) is identical to Problem (C9) which is merely a special case of Problem

PT

ED

(8). Hence, the assumption that (14) is true leads us to (C8) which confirms the existence of   Λ b, n ≤ v? < ∞ such that the solution to Problem (15) is the same as that of Problem (8).

AC

CE

Appendix D. Proof of Proposition 4 ( ) cn r Let p < 1 − knn where n = arg min . Then, by Proposition 3, (C3) and p > 0, we obtain n∈N rn −1 Λ(b, n) ≤ v? = FP

n∈N>0 (x? )e θn ( xn? )

−1 (p) where FP

n∈N>0 (x? )e θn ( xn? )

−1 Proposition 4 claims that v? = FP

n∈N>0 (x? )e θn ( xn? )

(p) > 0.

(D1)

(p) > 0 has a finite computable upper bound

† which for our purposes we shall call v[0] U B . Now, postulate the existence of a decision vector x

and a scalar α† both of which directly influence v[0] U B . In other words, suppose that there exists 41

ACCEPTED MANUSCRIPT

some scalar quantity   † ? −1 P v[0] U B x , α† ≥ v = F

n∈N>0 (x? )e θn ( xn? )

(p) > 0.

(D2)

                  

  † −1 v[0] (p) U B x , α† , F ψ e(x† )     X   h  i †   , e e x† = E  θ x µ† , E ψ n n    n∈N>0 (x† )

(D3)

AN US

and,

CR IP T

  e x† > 0 with, Moreover, associate to x† , a hypothetical random variable ψ

M

    X   h  i   2 †  e e x ≥ Var   σ† , Var ψ θn (xn ) for all x,       n∈N>0 (x)        =⇒ Fψe(x† ) (h (x)) ≤ FPn∈N (x) eθn (xn ) (h (x)) for all x, (by first order stochastic dominance) >0        [0] † −1 −1   =⇒ vU B x , α† = Fψe x† (p) ≥ FP (p) for all x,  e  ( )  n∈N>0 (x) θn (xn )        [0] −1   (p) = v? . =⇒ vU B x† , α† = Fψe−1(x† ) (p) ≥ FP e ? n∈N>0 (x? ) θn ( xn )

(D4)

ED

  Combining (D2), (D3) and (D4) with the fact that e θn xn† is an exponential random variable X e with rate parameter r 1x† , and θn (xn ) is the superposition of |N>0 (x)| uncorrelated expon n

n∈N>0 (x)

PT

nential random variables (each with rate parameter

1 rn xn ),

it becomes evident that a good way to

AC

CE

begin constructing the proof of Proposition 4 is to demonstrate the existence of an actual random   e x† > 0 such that variable ψ  X h  i    e x† = µ† = E ψ rn xn†       n∈N>0 (x† )   X h  i   2  e x† ≥  (rn xn )2 for all x σ = Var ψ    † n∈N>0 (x)

42

(D5)

ACCEPTED MANUSCRIPT

holds true. To show the veracity of (D5), we start by defining for all x

>0

>0

CR IP T

   X    X     =  e  (x ) rn xn θ µ , E  n n       n∈N>0 (x) n∈N >0 (x)     X X   X h i h i    2 2  =  e e e e  (x ) (r ) (x ) (x ) (x )  Var θ = x (since Cov θ , θ = 0 for all n , m) θ σ , Var n n n n n n m m n n     n∈N (x) n∈N (x) n∈N (x) >0

(D6)

and we note that a generalization of the Chebyshev lower bound on FPn∈N (x) eθn (xn ) (h (x)) = >0     X e θn (xn ) ≤ h (x) yields Pr  n∈N>0 (x)

AN US

2  σ2 + µ − h(x) 2 FPn∈N (x) eθn (xn ) (h (x)) ≥ 1 − .  h(x) 2 >0

(D7)

2

(

Next, introduce a constant α > 2 (that generalizes α† ), impose a new restriction σ2 = P 2 ) n∈N>0 (x) rn xn for all x, and then consider the following nonlinear optimization problem (α−1)(α−2)

1− α1

M

X     z x? ≤ y x† , α† , µ† , σ2† , min cn xn x

(D8)

AC

CE

PT

ED

α>2 n∈N µ,σ2

subject to  2 σ2 + µ − h(x) rn 2 1− ≥ p where p < 1 −  h(x) 2 kn

µ=

X

(D9)

2

rn xn (thanks to µ in (D6))

(D10)

n∈N>0 (x)

σ2 =

X

(rn xn )2 (thanks to σ2 in (D6))

(D11)

n∈N>0 (x)

σ2 =

 P 2  1 − α1 n∈N>0 (x) rn xn (α − 1)(α − 2)

xn ≥ ln for all n ∈ N.

(D12) (D13)

 
(D7) and (D9) both imply (9), and Problem (D8) is able to assert that z x? ≤ y x† , α† , µ† , σ2†

because of the additional restrictions (D10), (D11), (D12). By virtue of (D9), (D11), (D12) and 43

ACCEPTED MANUSCRIPT

the fact that Problem (D8) is a minimization problem, it follows that 

1−

1 α†

 P

 † 2 n∈N>0 (x† ) rn xn

>0

X β† µ† (rn xn )2 for all x where µ† =

≥ σ2 = ⇐⇒ σ2† = α† − 1 α† − 2 n∈N (x) >0

n∈N>0 (x† )

 rn xn† and β† = 1 −

α† β† > 0. α† − 1

Hence, σ2† = ⇐⇒ σ2† =

β†

α β  † †

α† −1




α† − 1 α† − 2 α† β2†
2

α† − 1

1 α†



X

rn xn† > 0

n∈N>0 (

x†

)

rn xn† (assuming (D5)).

n∈N>0 (x† )

(D14)

rn xn† > 0 both imply

AN US

µ† =

X

X

n∈N>0 (x† )


≥ σ2 =

X

(D15)

(rn xn )2 for all x, by (D14), (D15)

n∈N>0 (x)

X 2 (rn xn )2 for all x.
≥σ = α† − 2 n∈N>0 (x)

M

But µ† =

X

CR IP T

X (rn xn )2 for all x,

≥ σ2 = α† − 1 α† − 2 n∈N>0 (x) P ! X β† n∈N>0 (x† ) rn xn† 1 2 2 2 (rn xn ) for all x where β† , 1 −

≥σ = ⇐⇒ σ† = α† − 1 α† − 2 α† n∈N (x)

σ2† =

(D16)

ED

Recall that (D3) and (D4) respectively defined µ† and σ2† as the mean and variance of a   e x† . However, µ† in (D15) and σ2 in (D16) are respectively the hypothetical random variable ψ †

PT

mean and variance formulas for what is known as a Pareto random variable with shape parameter α† and scale parameter β† . In (D14), we used the expression for µ† in (D5) to derive (D15), and

AC

CE

note that (D16) merely conveys the same information as σ2† in (D5). We therefore conclude   e x† is an actual (Pareto) random variable satisfying the conditions laid out in (D5). The that ψ  P 1− α1 r x†   n∈N>0 (x† ) n n † β † e x† is F −1 † (p) = = and invoking (D4) we quantile function of ψ 1 1 e(x ) ψ (1−p) α† (1−p) α† obtain  1−   † −1 v[0] e(x† ) (p) = U B x , α† = F ψ

1 α†

P

† n∈N>0 (x† ) rn xn

(1 − p)

1 α†

−1 ≥ FP

e ( xn? ) (p)

n∈N>0 (x? ) θn

= v?

(D17)

where x† belongs to the optimal solution set for Problem (D8). Unfortunately, Problem (D8) is a difficult problem to solve mainly because of (D12). To surmount this obstacle, we need for 44

ACCEPTED MANUSCRIPT

ˆ , αˆ that satisfy v[0] U B (·) a pair of easily-computable arguments x

X

n∈N



=

1−

1 αˆ

P

ˆn n∈N>0 (ˆx) rn x

(1 − p)

  † ≥ v[0] U B x , α†

1 αˆ

where xˆ , αˆ need only be feasible for Problem (D8) since

(D18)

∂ [0] (ˆ ˆ ∂ xˆn vU B x, α)

=

1− α1ˆ

rn > 0 =⇒

CR IP T

v[0] x, α) ˆ U B (ˆ

1

(1−p) αˆ

cn xˆn ≥ y(x† , α† , µ† , σ2† ) for all n ∈ N. A direct way to determine xˆ , αˆ is to express them as

optimal solution set members for the following over-constrained version of Problem (D8): X       cn xn z x? ≤ y x† , α† , µ† , σ2† ≤ w xˆ , α, ˆ µ, ˆ σ ˆ 2 , min x

(D20)

AN US

α>2 n∈N µ,σ2

subject to  2 σ2 + µ − h(x) rn 2 1− = p where p < 1 − (D21)  h(x) 2 kn 2

M

h (x) µ= 2 X µ= rn xn

ED PT

(D23)

n∈N>0 (x)

σ2 =

X

(rn xn )2

n∈N>0 (x) 2

(D22)

σ =



1−

1 α

 P

n∈N>0 (x) rn xn

(α − 1)(α − 2)

xn ≥ ln for all n ∈ N.

(D24) 2

(D25) (D26)

CE

Then, using

h (ˆx) σ ˆ2 ⇐⇒ 2 = 1 − p, (D27) 2 µˆ   X 1 − α1ˆ µˆ 2 1 − α1ˆ σ ˆ2 2 (D23) and (D25): µˆ = rn xˆn ⇐⇒ σ ˆ = ⇐⇒ 2 = , (αˆ − 1)(αˆ − 2) (αˆ − 1)(αˆ − 2) µˆ n∈N (ˆx)

AC

(D21) and (D22): µˆ =

(D22) and (D23): µˆ =

>0

(D28)

X

h (ˆx) = 2 n∈N

x) >0 (ˆ

rn xˆn ⇐⇒ h (ˆx) = 2

45

X

n∈N>0 (ˆx)

rn xˆn

(D29)

ACCEPTED MANUSCRIPT

we obtain via (D27) and (D28), the quadratic equation 1 − α1ˆ 1 1 = ⇐⇒ αˆ 2 − 2αˆ − =0 (αˆ − 1)(αˆ − 2) αˆ (αˆ − 2) 1− p s 1 =⇒ αˆ = 1 + 1 + > 2 (by the quadratic formula). 1− p

(D30) (D31)

CR IP T

1− p=

whereas (6), Problem (D20), (D26) and (D29) give us the linear program X   cn xn z x? ≤ w (ˆx) , min x

(D32)

n∈N

subject to

! X rn rn xn log rn xn ≥2 kn (x) n∈N (x)

AN US

X

h (x) = log(1 − b) −

n∈N>0

(D33)

>0

xn ≥ ln for all n ∈ N.

(D34)

However, the optimal solution for Problem (D32) might be infeasible for Problem (D20). To

M

resolve this feasibility issue, we recall that v? is an ingredient of Problem (15) which is predicated on the validity of (14). But (14) is itself a generalization of (11), and Proposition 2 confirms that (12) is the optimal solution for Problem (8) when (11) is true. Consequently, to ensure that the

ED

optimal solution for Problem (D32) is feasible for Problem (D20), we embed (12) into Problem (D32) by replacing Problem (D32) with

PT

X     cn xn z x? ≤ y x† , α† , µ† , σ2† ≤ u (ˆx) , min

CE

x

(D35)

n∈N

subject to X h (x) ≥ 2 rn ln

(D36)

xn ≥ ln +

(D37)

AC

n∈N

log (1 − b)   for all n ∈ N. r 1 rn log knn · 1−p

Combining (D17), and (D18), we obtain

v? ≤



1−

1 αˆ

P

ˆn n∈N>0 (ˆx) rn x 1

(1 − p) αˆ 46

? = v[0] U B , vU B

(D38)

ACCEPTED MANUSCRIPT

where αˆ satisfying (D30) is given by (D31), and xˆ is the optimal solution for Problem (D35).

CR IP T

Appendix E. Proof of Proposition 5 Let s be a complex number with real part <(s) and for n ∈ N>0 (x), denote by feθn (xn ) (·) the

probability density function of an exponential random variable with rate parameter rn1xn . AddiR h(x) tionally, let Feθn (xn ) (h (x)) , 0 feθn (xn ) (t) dt indicate the cumulative distribution of e θn (xn ) evaluated at h (x) ≥ 0. The Laplace transform of feθn (xn ) (h (x)) is the following integral:

AN US

Z ∞ n o L feθn (xn ) (h (x)) (s) , exp (−sh (x)) · feθn (xn ) (h (x)) dh (x) 0 ! Z ∞ h (x) 1 exp − dh (x) = exp (−sh (x)) · rn xn rn xn 0 ( ) 1 1 = where <(s) > max − n∈N>0 (x) 1 + rn xn s rn xn

n o and h (x) ≥ 0. Since all the members in e θn (xn )

n∈N>0 (x)

M

implies

(E1)

are mutually-uncorrelated, (E1)

(E2)

AC

CE

PT

ED

  g(s) , L fPn∈N (x) eθn (xn ) (h (x)) (s) >0   = L feθ1 (h (x)) ⊕ . . . ⊕ feθ N (x) (h (x)) (s) | >0 | Y 1 = 1 + rn xn s n∈N>0 (x) Y Y 1 1 = · r x s + rn1xn n∈N>0 (x) n n n∈N>0 (x) Y Y 1 1 = ·  |Mi(n) (x)| r x 1 n∈N>0 (x) n n i(n)∈I(x) s + ri(n) xi(n)

where <(s) > max

i(n)∈I(x)

(

(E3) (E4)

) 1 − , and ⊕ stands for the convolution operation. Applying the ri(n) xi(n)

47

ACCEPTED MANUSCRIPT

residue method of partial fraction decomposition, we obtain g(s) 1 n∈N>0 (x) rn xn

=

Y

i(n)∈I(x)

=

1 

1 ri(n) xi(n)

i(n) (x)| X |MX

j=1

i(n)∈I(x)

=

s+

X

i(n)∈I(x)



∆i(n), j

s+

∆i(n),1 s+

|Mi(n) (x)| by (E4)

1 ri(n) xi(n)

j 1 ri(n) xi(n)

+

∆i(n),2

s+

where  ∆i(n), j =  Mi(n) (x) − j ! 1  =  Mi(n) (x) − j !

lim

d|Mi(n) (x)|− j ds|Mi(n) (x)|− j

lim

d|Mi(n) (x)|− j ds|Mi(n) (x)|− j

1 i(n) xi(n)

s→− r

1 i(n) xi(n)

s→− r

2 + . . . + 

∆i(n),|Mi(n) (x)| |Mi(n) (x)| s + ri(n)1xi(n)

(E5)

  !|Mi(n) (x)|   1 g(s)  s +  Q   1 ri(n) xi(n) n∈N>0 (x) rn xn   !|Mi(n) (x)| Y   1 1  (by (E3))  s + 1 ri(n) xi(n) s+ r x  n∈N (x)

AN US

1

1 ri(n) xi(n)

CR IP T

Q

>0

n n

(E6)

AC

CE

PT

ED

M

n o for i(n) ∈ I (x), and j ∈ 1, . . . , Mi(n) (x) . Referring back to (E2) and taking the inverse

48

ACCEPTED MANUSCRIPT

Laplace transform on both sides, we obtain for h (x) ≥ 0, fPn∈N

e

θ (x ) >0 (x) n n

(h (x))

=

i(n) (x)| X |MX

j=1

h j−1 (x) h (x) exp − ( r x j − 1)! r i(n) xi(n) n∈N>0 (x) n n ∆i(n), j

PT

i(n)∈I(x)

Q

ED

M

AN US

CR IP T

= L−1 {g(s)} (h (x))       Y   g(s) 1   (h (x)) = L−1  ·  Q  1    n∈N>0 (x) rn xn n∈N>0 (x) rn xn      Y   1 −1  g(s)   (h (x)) = L  Q   1    r x n∈N>0 (x) rn xn n∈N>0 (x) n n         Y X ∆   ∆i(n),1 ∆i(n),2 i(n),|Mi(n) (x)| 1 −1   (h (x)) by (E5) + . . . + L  + =        1 2 (x) M   | | i(n) r x 1   s + 1  i(n)∈I(x) s + ri(n) xi(n) ri(n) xi(n) n∈N>0 (x) n n s + ri(n) xi(n)   ! !  Y 1   X h (x) h (x) h (x)   =  ∆i(n),1 exp − + ∆i(n),2 exp − r x i(n)∈I(x) ri(n) xi(n) 1! ri(n) xi(n) n∈N>0 (x) n n  ! h|Mi(n) (x)|−1 (x) h (x)   exp − + . . . + ∆i(n),|Mi(n) (x)|   ri(n) xi(n)  Mi(n) (x) − 1 ! ! ! X ∆i(n),1 ∆i(n),2 h (x) h (x) h (x) Q exp − = +Q exp − ri(n) xi(n) ri(n) xi(n) n∈N>0 (x) rn xn n∈N>0 (x) rn xn 1! i(n)∈I(x) ! ∆i(n),|Mi(n) (x)| h|Mi(n) (x)|−1 (x) h (x)  exp −  + ... + Q ri(n) xi(n) n∈N>0 (x) rn xn Mi(n) (x) − 1 ! !

AC

CE

where ∆i(n), j is given by (E6). The cumulative probability measure FPn∈N

49

(E7)

e

θ (x ) >0 (x) n n

(h (x)) is

ACCEPTED MANUSCRIPT

given by the integral

e θ (x ) >0 (x) n n

(h (x)) =

Z

h(x)

0

=

Z

h(x)

0

=

fPn∈N

i(n) (x)| X |MX

j=1

i(n) (x)| X |MX

i(n)∈I(x)

(t) dt

i(n) (x)| X |MX

j=1

i(n)∈I(x)

i(n)∈I(x)

=

e

θ (x ) >0 (x) n n

j=1

Q

(E8)

! t t j−1 exp − dt by (E7) ri(n) xi(n) n∈N>0 (x) rn xn ( j − 1)! ∆i(n), j

1 n∈N>0 (x) rn xn ( j − 1)!

Q

∆i(n), j

Z

0

CR IP T

FPn∈N

h(x)

t j−1 exp −

! t dt ri(n) xi(n)

1 h j (x) h (x) Q  j Γ ( j) − Γ j,  h(x) ( r x j − 1)! r i(n) xi(n) n∈N>0 (x) n n ∆i(n), j

ri(n) xi(n)

AN US

   i(n) (x)|   X |MX Γ j, ri(n)h(x) ∆i(n), j x i(n)  1 −  r x
j Q = i(n) i(n)   ( j − 1)!  n∈N>0 (x) rn xn j=1 i(n)∈I(x)

!!

(E9)

  R∞ where Γ ( j) , ( j − 1)! = Γ ( j, 0), and Γ j, ri(n)h(x) t j−1 exp(−t)dt are the Gamma, h(x) xi(n) , ri(n) xi(n)

and upper incomplete Gamma function respectively. Finally, note that FPn∈N

e

θ (x ) >0 (x) n n

(h (x)) = 0

for h (x) < 0 since the random variable resulting from the superposition of a set of exponential

AC

CE

PT

ED

M

random variables can never take on a negative value.

50