DATA-DRIVEN DISTRIBUTIONALLY ROBUST OPTIMIZATION

Abstract

Embodiments of the disclosure include a method for providing data-driven distributionally robust optimization. The method includes receiving a plurality of samples of one or more uncertain parameters for a complex system and calculating a distribution uncertainty set for the one or more uncertain parameters. The method also includes receiving a deterministic problem model associated with the complex system that includes an objective and one or more constraints and creating a distributionally robust counterpart (DRC) model based on the distribution uncertainty set and the deterministic problem model. The method further includes formulating the DRC as a generalized problem of moments (GPM), applying a semi-definite programing (SDP) relaxation to the GPM and generating an approximation for a globally optimal distributionally robust solution to the complex system.

Description

CROSS-REFERENCE TO RELATED APPLICATION

The present application is a continuation of and claims priority from U.S. patent application Ser. No. 13/587,980 filed on Aug. 17, 2012 the entire contents of which are incorporated herein by reference.

BACKGROUND

The present invention relates to optimization of complex systems under uncertainty, and more specifically, to data-driven distributionally robust optimization of complex systems.

Many of today's complex systems require decision making that is affected by uncertainty in one or more system parameters, such as usage or demand. Currently, data relating to the uncertain aspects of these systems is periodically collected by one or more sensors or meters. Current robust optimization models only exploit the collected usage data for support information of distributions for the uncertain parameters, which often leads to overly conservative models. Current distributionally robust optimization systems only exploit the observed data to construct distributional uncertainty sets consistent with the first two moments of the observed data and/or handle only restrictive classes of objective functions and constraints. On the other hand, current stochastic optimization models require highly accurate knowledge of distribution of uncertain system parameters.

Furthermore, in large-scab: systems characterization of uncertainty of system parameters given based on collected data can be challenging.

Optimization models of complex real-world systems often involve nonlinear functionalities in objective and constraints. However, current distributionally robust optimization models can not take into account broad classes of nonlinearities in objective and constraints.

SUMMARY

Embodiments include a method for providing data-driven distributionally robust optimization. The method includes receiving a plurality of samples of one or more uncertain parameters for a complex system, the plurality of samples being described by an unknown distribution, and calculating a distribution uncertainty set for the one or more uncertain parameters. The method also includes receiving a deterministic problem model associated with the complex system that includes an objective and one or more constraints and creating a distributionally robust counterpart (DRC) model based on the distributional uncertainty set and the deterministic problem model. The method further includes formulating the DRC as a generalized problem of moments (GPM), applying a semi-definite programing (SDP) relaxation to the GPM and generating an approximation for a globally optimal distributionally robust solution to the complex system.

Additional features and advantages are realized through the techniques of the present invention. Other embodiments and aspects of the invention are described in detail herein and are considered a part of the claimed invention. For a better understanding of the invention with the advantages and the features, refer to the description and to the drawings.

BRIEF DESCRIPTION OF THE SEVERAL VIEWS OF THE DRAWINGS

The subject matter which is regarded as the invention is particularly pointed out and distinctly claimed in the claims at the conclusion of the specification. The forgoing and other features, and advantages of the invention are apparent from the following detailed description taken in conjunction with the accompanying drawings in which:

FIG. 1 illustrates a block diagram of a computer system for use in practicing the teachings herein;

FIG. 2 illustrates a flow diagram of a method for providing data-driven distributionally robust optimization in accordance with an embodiment; and

FIG. 3 illustrates a block diagram of a system for providing data-driven distributionally robust optimization in accordance with an embodiment.

DETAILED DESCRIPTION

In accordance with exemplary embodiments, an optimization model is provided that is robust with respect to uncertain input data, but not too conservative due to extreme scenarios, and is able to handle a variety of nonlinear objectives and constraints.

In accordance with exemplary embodiments, systems and computer program products for deriving approximations of globally optimal, distributionally robust solutions to optimization problems with parameter uncertainty in the presence of samples of the uncertain parameters are provided. In exemplary embodiments, the system takes into account uncertainty in the input data and employs convex optimization solvers to provide globally optimal distributionally robust solutions for a broad class of models for real-world systems.

FIG. 1 illustrates a block diagram of a computer system 100 for use in practicing the teachings herein. The methods described herein can be implemented in hardware, software (e.g., firmware), or a combination thereof. In an exemplary embodiment, the methods described herein are implemented in hardware, and may be part of the microprocessor of a special or general-purpose digital computer, such as a personal computer, workstation, minicomputer, or mainframe computer. The computer system 100 therefore includes general-purpose computer 101.

In an exemplary embodiment, in terms of hardware architecture, as shown in FIG. 1, the computer 101 includes a processor 105, memory 110 coupled to a memory controller 115, and one or more input and/or output (I/O) devices 140, 145 (or peripherals) that are communicatively coupled via a local input/output controller 135. The input/output controller 135 can be, for example but not limited to, one or more buses or other wired or wireless connections, as is known in the art. The input/output controller 135 may have additional elements, which are omitted for simplicity, such as controllers, buffers (caches), drivers, repeaters, and receivers, to enable communications. Further, the local interface may include address, control, and/or data connections to enable appropriate communications among the aforementioned components.

The processor 105 is a hardware device for executing hardware instructions or software, particularly that stored in memory 110. The processor 105 can be any custom made or commercially available processor, a central processing unit (CPU), an auxiliary processor among several processors associated with the computer 101, a semiconductor based microprocessor (in the form of a microchip or chip set), a macroprocessor, or generally any device for executing instructions. The processor 105 includes a cache 170, which may include, but is not limited to, an instruction cache to speed up executable instruction fetch, a data cache to speed up data fetch and store, and a translation lookaside buffer (TLB) used to speed up virtual-to-physical address translation for both executable instructions and data. The cache 170 may be organized as a hierarchy of more cache levels (L1, L2, etc.).

The memory 110 can include any one or combination of volatile memory elements (e.g., random access memory (RAM, such as DRAM, SRAM, SDRAM, etc.)) and nonvolatile memory elements (e.g., ROM, erasable programmable read only memory (EPROM), electronically erasable programmable read only memory (EEPROM), programmable read only memory (PROM), tape, compact disc read only memory (CD-ROM), disk, diskette, cartridge, cassette or the like, etc.). Moreover, the memory 110 may incorporate electronic, magnetic, optical, and/or other types of storage media. Note that the memory 110 can have a distributed architecture, where various components are situated remote from one another, but can be accessed by the processor 105.

The instructions in memory 110 may include one or more separate programs, each of which comprises an ordered listing of executable instructions for implementing logical functions. In the example of FIG. 1, the instructions in the memory 110 include a suitable operating system (OS) 111. The operating system 111 essentially controls the execution of other computer programs and provides scheduling, input-output control, file and data management, memory management, and communication control and related services.

In an exemplary embodiment, a conventional keyboard 150 and mouse 155 can be coupled to the input/output controller 135. Other output devices such as the I/O devices 140, 145 may include input devices, for example but not limited to a printer, a scanner, microphone, and the like. Finally, the I/O devices 140, 145 may further include devices that communicate both inputs and outputs, for instance but not limited to, a network interface card (NIC) or modulator/demodulator (for accessing other files, devices, systems, or a network), a radio frequency (RF) or other transceiver, a telephonic interface, a bridge, a router, and the like. The system 100 can further include a display controller 125 coupled to a display 130. In an exemplary embodiment, the system 100 can further include a network interface 160 for coupling to a network 165. The network 165 can be an IP-based network for communication between the computer 101 and any external server, client and the like via a broadband connection. The network 165 transmits and receives data between the computer 101 and external systems. In an exemplary embodiment, network 165 can be a managed IP network administered by a service provider. The network 165 may be implemented in a wireless fashion, e.g., using wireless protocols and technologies, such as WiFi, WiMax, etc. The network 165 can also be a packet-switched network such as a local area network, wide area network, metropolitan area network, Internet network, or other similar type of network environment. The network 165 may be a fixed wireless network, a wireless local area network (LAN), a wireless wide area network (WAN) a personal area network (PAN), a virtual private network (VPN), intranet or other suitable network system and includes equipment for receiving and transmitting signals.

If the computer 101 is a PC, workstation, intelligent device or the like, the instructions in the memory 110 may further include a basic input output system (BIOS) (omitted for simplicity). The BIOS is a set of essential routines that initialize and test hardware at startup, start the OS 111, and support the transfer of data among the hardware devices. The BIOS is stored in ROM so that the BIOS can be executed when the computer 101 is activated. When the computer 101 is in operation, the processor 105 is configured to execute instructions stored within the memory 110, to communicate data to and from the memory 110, and to generally control operations of the computer 101 pursuant to the instructions.

In exemplary embodiments, a generalized problem of moments (GPM) with polynomial data is an optimization problem of the form:

mi_+(K)∫p(x)dμ(x)

s.t. ∫f_i(x)dμ(x)=b_i, iεl,

∫f_j(x)dμ(x)≦b_j, jεJ,   (1)

where K={xεⁿ|g₁(x)≧0, . . . , g_m(x)≧0}, p, f_i(iεl), f_j(jεJ)_gj(j=1, . . . , m)ε[x], I, J a finite or infinite set of indices, and ₊(K)a set of finite Borel measure supported on K. The GPM is an infinitely dimensional linear program. It can be approximated by the Lasserre's hierarchy of SDP relaxations:

min_y L_y(p)

s.t. M_d(y)0,

M_d−dj (g_jy)0, j=1, . . . , m,

L_y(f_i)=b_i, iεl,

L_y(f_i)≦b_j, jεJ,

y₀≦Const.,   (2)

where

$d_j:=\lceil \frac{\mathrm{deg}\left(g_j\right)}{2}\rceil$

and dεN the order of a single SDP relaxation. Let min (GPM) the minimum of (1), min (SDP_d) the minimum and y_d* the minimizer of (2). If K compact and its quadratic module archimedian, min (SDP_d)→min (GPM) for d→∞. Moreover, if (1) has a unique measure as minimal solution, then y_d* converges to the moment vector of the optimal measure for (1) as well. In exemplary embodiments, M_d(y) and M_d−dj(g_ly) are linear combinations of the components y_αof y with real-valued, symmetric matrices as multipliers, and M_d(y) stands for a constraint where the matrix M_d(y) is required to be positive semidefinite. In exemplary embodiments L_y(p) is a linear operator that maps a polynomial P(x)=E_αP_αx_α to a linear combinations of the components y_α of y such that L_y(p)=Σ_αP_αyα. Note, if μεM₊(K) is additionally restricted to be a probability measure, the GPM (1) becomes a Polynomial Optimization Problem (POP). Thus, in exemplary embodiments, a POP can be approximated by Lasserre's hierarchy of SDP (2) as well.

In exemplary embodiments when no knowledge of a distribution of a uncertain parameter over the uncertainty set is assumed, the optimization problem can be represented as:

min_xεX maE_μ(ξ)[h(x,ξ)]=min_xεX ma∫_Sh(x,ε)dμ(ξ)=min_xεX maF(x,μ)   (3)

where some set of finite Borel measures supported on S⊂^m, as the distributionally robust counterpart (DRC) of the problem:

min_xεX h(x, ξ)  (4)

for some given ξεS. In the data-driven case, is constructed from a given sample ξ⁽¹⁾, . . . , ξ^(N) of the actual, unknown distribution μ* of the uncertain parameter ξ.

In exemplary embodiments, an uncertainty set is calculated from a given sample of collected data. In exemplary embodiments, uncertainty sets can be defined by small deviations around the statistical moments. The statistical moments of μ* are given by:

$m_{\alpha }^N=\frac{1}{N}\sum _{i=1}^N{\left({\xi }^{\left(i\right)}\right)}^{\alpha },$

for all αε^m. The uncertainty set _d,ε,N can be represented as:

_d,ε,N={με₊(S)|∫_Sξ^αdμ(ξ)≦m_α^N+ε_α, ∫_Sξ^αdμ(ξ)≧m_α^N−ε_α∀|α|≦d},

for some dε and ε_α>0, where ₊(S) the set of finite Borel measures supported on the set S. ε needs to be chosen such that μ* ε. Furthermore, assuming the support S:={ξ|g_j(ξ)≧0(j=1, . . . , r)} of μ* where g_j ε[ξ] are multivariate real polynomials, is compact, and hε[x,ξ], i.e. a multivariate polynomial in x and ξ.

In exemplary embodiments, the distributionally robust counterpart model can be reformulated as generalized problem of moments. Utilizing the dual of the inner maximization problem, (3) can be re-written as:

min_xεX,λ⁺≦0,λ⁻≦0 Σ_αλ_α⁺b_α⁺+Σ_aλ_α⁻b_α⁻

s.t. Σ_αλ_α⁺h_α⁺(ξ)+Σ_αλ_α⁻h_α⁻(ξ)≦h(x, ξ)∀ξεS,   (5)

where h_α⁺(ξ):=−ξ^α, h_α⁻(ξ):=ξ^α, b_α⁺:=−m_α^N−ε_{α and bα}⁻:=m_α^N−ε_α, i.e. ∫_Sh_α⁺(ξ)dμ(ξ)≧b_α⁺and ∫_Sh_α⁻(ξ)dμ(ξ)≧b_α⁻ the inequalities in the definition of _d,ε,N. Assuming h(x,ξ):=Σ_αh_α(ξ)x^α where h_αε[ξ]. Define z:=(x,λ⁺, λ⁻) and Z:=X×₋^2k. Then, (5) can be rewritten as

$\begin{array}{cc}\underset{z\in Z}{\min }\sum _{\alpha }{\tilde{b}}_{\alpha }z^{\alpha }s.t.\sum _{\beta }{\tilde{h}}_{\beta }\left(\xi \right)z^{\beta }\ge 0\forall \xi \in S,& \left(6\right)\end{array}$

where {tilde over (b)} and {tilde over (h)}_β defined based on b, h_α⁺, h_α⁻ and h. Then, the polynomial constraint (6) can be tightened by introducing r matrix variables A_j:

$\begin{array}{cc}\underset{z\in Z,A_j\in {}^{s\left(d\right)}}{\min }\sum _{\alpha }{\tilde{b}}_{\alpha }z^{\alpha }s.t.\begin{array}{cc}\sum _{\beta }{\tilde{h}}_{\beta }\left(\xi \right)z^{\beta }=\sum _{j=0}^rg_j\left(\xi \right)u_d\left(\xi \right)A_ju_d\left(\xi \right)& \forall \xi \in {ℝ}^m,\\ A_j\succcurlyeq 0& \forall j\in \left\{0,\dots ,r\right\},\end{array}& \left(7\right)\end{array}$

where g₀(ξ):=1 and u_d(ξ)=(1,ξ₁, . . . , ξ_m^d) the standard basis of [ξ]_d.

The set of polynomial equality constraints in (7) can be rewritten as a set of constraints linear in z^β and components of A_j. Since each of the constraints A_jε₊^s(d) can be written as a finite number of scalar, polynomial inequality constraints, (7) is equivalent to a polynomial optimization problem, which can be approximated by a sequence of semi-definite programing (SDP) relaxations.

In exemplary embodiments, the of (3) is a set of Borel measures defined by their densities with respect to the Lebesgue measure supported on the same set S as before. The statistical estimate for the density of μ* based on the sample ξ⁽¹⁾, . . . , ξ^(N) can be denote by f^NεL²(S). Any one of numerous known methods including, but not limited to, Kernel Density Estimation or wavelets, can be used to estimate densities of an unknown measure given a finite sample. In exemplary embodiments, f^N is a multivariate polynomial, in order to guarantee the set D_d,ε,N constructed below is nonempty when d is fixed.

In an exemplary embodiment, the uncertainty set may be represented as:

_ε,N={με₊(S)|∫_Sdμ(ξ)=∫_Sf(ξ)dξ,∫_S(f(ξ)−f^N(ξ))²dξ≦ε}.

In order to derive tractable distributionally robust counterpart model (DRC), we consider the following truncated, polynomial approximation for D_ε,N:

_d,ε,N={fε[ξ]_d|f(ξ)≧0∀ξεS,∫_S(f(ξ)−f^N(ξ))²dξ≦ε},   (8)

i.e. the inner maximization problem is equivalent to:

ma_[ξ]d Σ_αx^α∫_Sh_α(ξ)f(ξ)dξ

s.t. f(ξ)≧0 ∀ξεS,

∫_S(f(ξ)−f^N(ξ))²dξ≦ε.   (9)

For a fixed xεX, (9) is a polynomial optimization problem. Its dual is a minimization problem involving polynomials, moment expressions and the closed, semialgebraic set S as well. Therefore, the DRC (3) can be reformulated as a minimization problem involving moments and polynomials, which can be approximated by a converging sequence of SDP relaxations.

In another exemplary embodiment, the uncertainty set may be represented as:

_ε,N={με₊₍s)|∫_Sdμ(ξ)=∫_Sf(ξ)dξ, max_ξεS|f(ξ)−f^N(ξ)|≦ε}.

In order to derive tractable DRC, _ε,N can be approximated as:

_d,ε,N={fε[ξ]_d|f(ξ)≧0∀ξεS, max_ξεS|f(ξ)−f^N(ξ)|≦ε},   10)

i.e. the inner maximization problem is equivalent to:

ma_[ξ]d Σ_αx^α∫_Sh_α(ξ)f(ξ)dξ

s.t. f(ξ)≧0 ∀ξεS,

−f(ξ)+f^N(ξ)−ε≧0 ∀ξεS,

f(ξ)−f^N(ξ)+ε≧0 ∀ξεS.   (11)

With f(ξ)=Σ_βf_βξ^ξ, (11) can be rewritten as:

$\begin{array}{cc}\begin{array}{cccc}{\max }_{{\left(f_{\beta }\right)}_{{\beta }^{\in {ℝ}^q}}}& \sum _{\beta }^{}f_{\beta }h_{\beta }\left(x\right)-\sum _{\beta }^{}f_{\beta }{\xi }^{\beta }\le 0& \forall \xi \in S,& \\ s.t.& \sum _{\beta }^{}f_{\beta }{\xi }^{\beta }\le f^N\left(\xi \right)-\epsilon & \forall \xi \in S,& \\ & -\sum _{\beta }^{}f_{\beta }{\xi }^{\beta }\le -f^N\left(\xi \right)+\epsilon & \forall \xi \in S,& \end{array}& \left(12\right)\end{array}$

where h_β(x):=Σ_αx^αh_α,β:=Σ_αx^α∫_Sh_α(ξ)ξ^{βdξ. Taking the dual of (}12), we are able to reformulate the min-max problem (3) as the minimization problem:

$\begin{array}{cc}\begin{array}{ccc}{\min }_{x\in X,{\mu }_1\in {}_{+}\left(S\right),{\mu }_2\in {}_{+}\left(S\right),{\mu }_3\in {}_{+}\left(S\right)}& {\int }_S^{}0{\mu }_1\left(\xi \right)+{\int }_S^{}f^N\left(\xi \right)-\epsilon {\mu }_2\left(\xi \right)+{\int }_S^{}-f^N\left(\xi \right)+\epsilon {\mu }_3\left(\xi \right)& \\ s.t.& -{\int }_S^{}{\xi }^{\beta }{\mu }_1\left(\xi \right)+{\int }_S^{}{\xi }^{\beta }{\mu }_2\left(\xi \right)-{\int }_S^{}{\xi }^{\beta }{\mu }_3\left(\xi \right)=h_{\beta }\left(x\right)& \forall \beta .\end{array}& \left(13\right)\end{array}$

Assuming X:={x|k_j(x)≧0, j=1, . . . , t} compact, where k_jε[x]. (13) is equivalent to:

$\begin{array}{cc}\begin{array}{ccc}\underset{v\in {}_{+}\left(X\right),{\mu }_1\in {}_{+}\left(S\right),{\mu }_2\in {}_{+}\left(S\right),{\mu }_3\in {}_{+}\left(S\right)}{\min }& {\int }_S^{}0{\mu }_1\left(\xi \right)+{\int }_S^{}f^N\left(\xi \right)-\epsilon {\mu }_2\left(\xi \right)+{\int }_S^{}-f^N\left(\xi \right)+\epsilon {\mu }_3\left(\xi \right)& \\ s.t.& -{\int }_S^{}{\xi }^{\beta }{\mu }_1\left(\xi \right)+{\int }_S^{}{\xi }^{\beta }{\mu }_2\left(\xi \right)-{\int }_S^{}{\xi }^{\beta }{\mu }_3\left(\xi \right)={\int }_X^{}h_{\beta }\left(x\right)v\left(x\right)& \forall \beta ,\\ & {\int }_X^{}v\left(x\right)=1.& \end{array}& \left(14\right)\end{array}$

Problem (14) is a Generalized Problem of Moments (GPM) with polynomial data, whose minimum can be approximated up to arbitrary precision by a hierarchy of SDP relaxations, if X and S are compact and archimedian. Moreover, in the case (14) has a unique minimizer, the sequence of optimal solutions of the hierarchy of SDPs converges to this minimizer for increasing relaxation order.

In exemplary embodiments, it can be assumed that samples ξ⁽¹⁾, . . . , ξ^(N) for the uncertain parameter ξ, which takes values in a given bounded interval [A, B]⊂, which can be broken into K-intervals:

u₀, . . . , u_K−1,

such that:

|u_k|=|B−A|/K for all k=0, . . . , K−1.

Let m₀, . . . , m_K−1 denote the midpoints of the respective intervals, the empirical distribution {circumflex over (F)}_N,K can be defined as:

${\hat{F}}_{N,K}\left(k\right)=\frac{1}{N}\sum _{i=1}^N1_{\left[{\xi }^{\left(i\right)}\in u_k\right]}$ $\mathrm{for}\mathrm{all}k=0,\dots ,K-1.$

In exemplary embodiments, the uncertainty set may be approximated by the following optimization problem:

min_xεX max_γεUε∫_ξεSh(x,ξ)γ(ξ)dξ,   (15)

where:

$U_{\varepsilon }=\left\{\gamma \text{:}\left[A,B\right]\to {ℝ}_{+}\text{:}\int \gamma \left(\xi \right)\xi =1,\gamma \left(z\right)-\sum _{k=1}^K{\hat{F}}_{N,K}\left(k\right)1_{\left[z\in u_k\right]}\le \varepsilon \mathrm{for}\mathrm{all}z\right\}.$

Let d=(d₁, . . . , d_K) denote a vector in ^K. The optimization problem (15) can be approximated by the following:

min_xεX max_dεWεΣ_k=0^K−1h(x,m_k)d_k,   (16)

where the uncertainty set is:

$W_{\varepsilon }=\left\{d\in {\Delta }^K\text{:}d_k-{\hat{F}}_{N,K}\left(k\right)\le \varepsilon \frac{B-A}{K}\mathrm{for}\mathrm{all}k\right\}$

In exemplary embodiments, the optimization problem can be reformulated as a polynomial optimization problem. Observe that the inner maximization of (16) can be written as:

$\begin{array}{cc}{\max }_{d\in {ℝ}^K}\sum _{k=0}^{K-1}h\left(x,m_k\right)d_ks.t.d_k-{\hat{F}}_{N,K}\left(k\right)\le \varepsilon ,k=0,\dots ,K-1;{\hat{F}}_{N,K}\left(k\right)-d_k\le \varepsilon ,k=0,\dots ,K-1;\sum _{i=1}^Kd_i\le 1;d_j\ge 0,j=0,\dots ,K-1.& \left(M1\right)\end{array}$

The dual of the maximization (M1) is:

$\underset{y\in {ℝ}^{2K+1}}{\min }\sum _{k=0}^{K-1}y_k\left(\varepsilon +{\hat{F}}_{N,K}\left(k\right)\right)+\sum _{=K}^{2K-1}y_{}\left(\varepsilon -{\hat{F}}_{N,K}\left(-K\right)\right)+y_{2K}$ $s.t.y_k-y_{K+k}+y_{2K}\ge h\left(x,m_k\right),k=0,\dots ,K-1,y_j\ge 0,j=0,\dots ,2K.$

By the Duality Theorem for linear programs, the primal and dual have the same optimal value; hence, the optimization (16) can be written as:

$\begin{array}{cc}\underset{x\in X,y\in {ℝ}^{2K+1}}{\min }\sum _{k=0}^{K-1}y_k\left(\varepsilon +{\hat{F}}_{N,K}\left(k\right)\right)+\sum _{=K}^{2K-1}y_{}\left(\varepsilon -{\hat{F}}_{N,K}\left(-K\right)\right)+y_{2K}s.t.y_k-y_{K+k}+y_{2K}\ge h\left(x,m_k\right),k=0,\dots ,K-1,y_j\ge 0,j=0,\dots ,2K.& \left(A1\right)\end{array}$

(A1) is a polynomial optimization problem (POP) of dimension n+2K+1, its degree coincides with the degree of h.

If h is a multivariate polynomial in the arguments x and ξ, then the optimization problem:

min_xεX max_dεWεΣ_k=0^K−1h(x,m_k)d_k,   (17)

can be approximated by a sequence of SDP relaxations.

In exemplary embodiments, ξ denotes the uncertain parameter, which is a random variable takes values in a set S[A,B]^d⊂^d and ξ⁽¹⁾, . . . , ξ^(N) denotes a sequence of random variables with the same probability distribution as ξ. The set can be partitioned [A,B]^d into a regular grid of 10 hypercubes of equal volume:

{v_α:αε[K]^d}

Letting m_a denote the center of the hypercube v_a for every a, the partition of [A,B] can be defined as K intervals of equal length:

u₀, . . . , u_K−1.

ξ_l⁽ⁱ⁾ denotes the l-th component of the sample ξ⁽ⁱ⁾, for every i, and every pair l≠l′, the random variables ξ_l⁽ⁱ⁾ and ξ_l′⁽ⁱ⁾ are independent. If N and K are fixed; the marginal empirical frequencies {circumflex over (F)}_l can be defined as:

${\hat{F}}_{}\left(k\right)=\frac{1}{N}\sum _{i=1}^N1_{\left[{\xi }_{}^{\left(i\right)}\in u_k\right]}\mathrm{for}\mathrm{all}k=0,\dots ,K-1,\mathrm{and}$ $=1,\dots ,d.$

Let a₁, . . . , a_d denote the components of aε[K]^d. The joint empirical frequencies can be defined as:

$\hat{G}\left(a\right)=\prod _{=1}^d{\hat{F}}_{}\left(a_{}\right)\mathrm{for}\mathrm{all}a\in {\left[K\right]}^d.$

Consider the robust optimization problem with an uncertainty set centered on the empirical density:

min_xεX max_γεUε∫_ξεSh(x,ξ)γ(ξ)dξ,   (18)

where:

$U_{\varepsilon }=\left\{\gamma \text{:}S\to {ℝ}_{+}\text{:}\int \gamma \left(z\right)z=1,\gamma \left(z\right)-\sum _{a\in {\left[K\right]}^d}\hat{G}\left(a\right)1_{\left[z\in v_a\right]}\le \varepsilon \mathrm{for}\mathrm{all}z\right\}.$

The optimization problem (18) can be approximated by the following:

min_xεX max_pεWεΣ_αε[K]^dh(x,m_a)p(a),   (19)

where the uncertainty set is:

W_ε={p:[K]^d→₊:Σ_αp(α)≦1,|p(α)−Ĝ(α)|≦ε for all α}.

In exemplary embodiments, the optimization problem can be reformulated as a polynomial optimization problem. The inner maximization of (19) can be written as:

$\begin{array}{cc}\underset{x}{\min }\underset{p}{\max }\sum _{a\in {\left[K\right]}^d}^{}h\left(x,m_a\right)p\left(a\right)s.t.p\left(a\right)-\hat{G}\left(a\right)\le \varepsilon ,\forall a\in {\left[K\right]}^d;\hat{G}\left(a\right)-p\left(a\right)\le \varepsilon ,\forall a\in {\left[K\right]}^d;\sum _{a\in {\left[K\right]}^d}^{}p\left(a\right)\le 1;p\left(a\right)\ge 0,\forall a\in {\left[K\right]}^d.& \left(M2\right)\end{array}$

The dual of the maximization (M2) is:

$\underset{x}{\min }\underset{y,y^{\prime },y^{″}}{\min }\sum _{a\in {\left[K\right]}^d}^{}y_a\left(\varepsilon +\hat{G}\left(a\right)\right)+\sum _{a^{\prime }\in {\left[K\right]}^d}^{}y_{a^{\prime }}^{\prime }\left(\varepsilon -\hat{G}\left(a^{\prime }\right)\right)+y^{″}$ $s.t.y_a-y_a^{\prime }+y^{″}\ge h\left(x,m_a\right),\forall a\in {\left[K\right]}^d,y_a\ge 0,y_a^{\prime }\ge 0,y^{″}\ge 0,\forall a\in {\left[K\right]}^d.$

By the Duality Theorem for linear programs, the primal and dual have the same optimal value.

Referring now to FIG. 2, a flow diagram illustrating a method 200 for providing data-driven distributionally robust optimization in accordance with an embodiment is shown. As shown at block 202, the method 200 includes receiving sample of uncertain parameters for a system. In exemplary embodiments, the samples may be received from one or more sensors or meters in the system. Next, as shown at block 204, the method includes calculating a distribution uncertainty set using statistical tools based on polynomial probability density functions, moments up to arbitrary order or histogram estimates. In exemplary embodiments, the distributional uncertainty set may be described by polynomial inequality constraints. As shown at block 206, the method 200 also includes receiving a deterministic problem model that includes one objective and one or more constraints. In exemplary embodiments, the objective and the equality and inequality constraints may be described as multivariate real-valued polynomials. As shown at block 208, the method 200 includes creating a distributionally robust counterpart (DRC) model based on the distributional uncertainty set and the deterministic problem model. Next, as shown at block 210, the method 200 includes formulating the DRC as a Generalized Problem of Moments (GPM) with polynomial data. As shown at block 212, a Semi-Definite Programming (SDP) relaxation is applied to the GPM. The method 200 concludes at block 214 by generating an approximation for distributionally robust solution.

In exemplary embodiments, the approximation for distributionally robust solution includes a precision level that can be evaluated after the approximation for distributionally robust solution is created. In exemplary embodiments, the distributional uncertainty set includes a support that is described by one or more multivariate polynomial inequality constraints.

Referring now to FIG. 3, a block diagram of a system 300 for providing data-driven distributionally robust optimization is shown. As illustrated the system 300 includes a complex system 310 and an optimization system 304, which may be a computer system similar to the one shown in FIG. 1. The complex system 310 includes one or more sensors 312 that are configured to monitor the operation of one or more parameters of the complex system 310. In addition, the complex system 310 includes one or more controls 314 that are configured to control one or more operational characteristic of the complex system. In exemplary embodiments, the optimization system 304 receives data samples from the one or more sensors 312 and calculates an approximation for distributionally robust solution. Based on the approximation for distributionally robust solution, the controls 314 of the complex system 310 are configured to optimize the operation of the complex system 310.

In an exemplary embodiment, the complex system 310 may be a water distribution network characterized by a connected graph G (N, E), representing a water distribution network. Where Nis the set of nodes and E the set of pipes connecting the nodes. In addition, p_i the pressure, e_i elevation, and d_i demand at iεN, q_i,j the flow from i to j, and hl_i,j the headloss caused by friction in case of flow i to j for (i,j)εE. The optimization goal is to minimize the overall pressure in the water distribution network while adhering to the mass and energy conservations laws for flow and pressure:

min Σ_iεNp_i+Σ_jεN(d_j−Σ_k≠jq_k,j+Σ_l≠jq_j,l)²

s.t. p_min≦p_i≦p_max ∀iεN,

q_min≦q_i,j≦p_max ∀(i,j)εE,

q_i,j(p_j+e_j−p_i−e_i+hl_i,j)q_i,j))≦0 ∀(i,j)εE,

p_j+e_j−p_i−e_i+hl_i,j)q_i,j))≦0 ∀(i,j)εE,   (20)

Assuming the headloss, hl_i,j is a quadratic function in q and the vector of demands d=(d₁, . . . , d_|N|) is not known exactly and is therefore affected by uncertainty. Given a sample d⁽¹⁾, . . . , d^(K) of measurements of the demands at K discrete, equidistant time points.

Letting p=(p₁, . . . , p|N|)ε^{|N| and q=(q}_1,2, q_2,1, . . . ,)ε^{2|E| and defining h(p,q,d):=Σ}_iεNp_i+Σ_jεN(d_j−Σ_k≠jq_k,j+Σ_l≠jq_j,l)², X:={(p,q)ε^|N|+2|R||p_min≦p_i≦p_max ∀i, q_min≦q_i,j≦q_max, q_i,j(p_j+e_j−p_i−e_i+hl_i,j)≦0, p_j+e_j−p_i−e_i+hl_i,j(q_i,j)≧0∀(i,j)}. (20) is then equivalent to:

min_(p,q)εXh(p,q,d)   (21)

for a given demand profile d, i.e. it is of the form (4). In addition to the sampled demand, lower and upper bounds for demand at each node are given. Therefore,

dεS:={{tilde over (d)}ε^|N||d_i^min≦{tilde over (d)}_i≦d_i^max}.

The uncertainty set can be defined as:

_t,ε,K={fε[d]_t|f(d)≧0∀dεS, max_dεS|f(d)−f^K(d)|≦ε}.   (22)

The distributionally robust counterpart of (20) can be represented as:

$\underset{\left(p,q\right)\in X}{\min }\underset{f\in {}_{t,\epsilon ,K}}{\max }{\int }_Sh\left(p,q,d\right)f\left(d\right)d,$

i.e. it falls into the class (3).

As will be appreciated by one skilled in the art, aspects of the present invention may be embodied as a system, method or computer program product. Accordingly, aspects of the present invention may take the form of an entirely hardware embodiment, an entirely software embodiment (including firmware, resident software, micro-code, etc.) or an embodiment combining software and hardware aspects that may all generally be referred to herein as a “circuit,” “module” or “system.” Furthermore, aspects of the present invention may take the form of a computer program product embodied in one or more computer readable medium(s) having computer readable program code embodied thereon.

Any combination of one or more computer readable medium(s) may be utilized. The computer readable medium may be a computer readable signal medium or a computer readable storage medium. A computer readable storage medium may be, for example, but not limited to, an electronic, magnetic, optical, electromagnetic, infrared, or semiconductor system, apparatus, or device, or any suitable combination of the foregoing. More specific examples (a non-exhaustive list) of the computer readable storage medium would include the following: an electrical connection having one or more wires, a portable computer diskette, a hard disk, a random access memory (RAM), a read-only memory (ROM), an erasable programmable read-only memory (EPROM or Flash memory), an optical fiber, a portable compact disc read-only memory (CD-ROM), an optical storage device, a magnetic storage device, or any suitable combination of the foregoing. In the context of this document, a computer readable storage medium may be any tangible medium that can contain, or store a program for use by or in connection with an instruction execution system, apparatus, or device.

A computer readable signal medium may include a propagated data signal with computer readable program code embodied therein, for example, in baseband or as part of a carrier wave. Such a propagated signal may take any of a variety of forms, including, but not limited to, electro-magnetic, optical, or any suitable combination thereof. A computer readable signal medium may be any computer readable medium that is not a computer readable storage medium and that can communicate, propagate, or transport a program for use by or in connection with an instruction execution system, apparatus, or device.

Program code embodied on a computer readable medium may be transmitted using any appropriate medium, including but not limited to wireless, wireline, optical fiber cable, RF, etc., or any suitable combination of the foregoing.

Computer program code for carrying out operations for aspects of the present invention may be written in any combination of one or more programming languages, including an object oriented programming language such as Java, Smalltalk, C++ or the like and conventional procedural programming languages, such as the “C” programming language or similar programming languages. The program code may execute entirely on the user's computer, partly on the user's computer, as a stand-alone software package, partly on the user's computer and partly on a remote computer or entirely on the remote computer or server. In the latter scenario, the remote computer may be connected to the user's computer through any type of network, including a local area network (LAN) or a wide area network (WAN), or the connection may be made to an external computer (for example, through the Internet using an Internet Service Provider).

Aspects of the present invention are described below with reference to flowchart illustrations and/or block diagrams of methods, apparatus (systems) and computer program products according to embodiments of the invention. It will be understood that each block of the flowchart illustrations and/or block diagrams, and combinations of blocks in the flowchart illustrations and/or block diagrams, can be implemented by computer program instructions. These computer program instructions may be provided to a processor of a general purpose computer, special purpose computer, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, create means for implementing the functions/acts specified in the flowchart and/or block diagram block or blocks.

These computer program instructions may also be stored in a computer readable medium that can direct a computer, other programmable data processing apparatus, or other devices to function in a particular manner, such that the instructions stored in the computer readable medium produce an article of manufacture including instructions which implement the function/act specified in the flowchart and/or block diagram block or blocks.

The computer program instructions may also be loaded onto a computer, other programmable data processing apparatus, or other devices to cause a series of operational steps to be performed on the computer, other programmable apparatus or other devices to produce a computer implemented process such that the instructions which execute on the computer or other programmable apparatus provide processes for implementing the functions/acts specified in the flowchart and/or block diagram block or blocks.

The flowchart and block diagrams in the Figures illustrate the architecture, functionality, and operation of possible implementations of systems, methods and computer program products according to various embodiments of the present invention. In this regard, each block in the flowchart or block diagrams may represent a module, segment, or portion of code, which comprises one or more executable instructions for implementing the specified logical function(s). It should also be noted that, in some alternative implementations, the functions noted in the block may occur out of the order noted in the figures. For example, two blocks shown in succession may, in fact, be executed substantially concurrently, or the blocks may sometimes be executed in the reverse order, depending upon the functionality involved. It will also be noted that each block of the block diagrams and/or flowchart illustration, and combinations of blocks in the block diagrams and/or flowchart illustration, can be implemented by special purpose hardware-based systems that perform the specified functions or acts, or combinations of special purpose hardware and computer instructions.

The terminology used herein is for the purpose of describing particular embodiments only and is not intended to be limiting of the invention. As used herein, the singular forms “a”, “an” and “the” are intended to include the plural forms as well, unless the context clearly indicates otherwise. It will be further understood that the terms “comprises” and/or “comprising,” when used in this specification, specify the presence of stated features, integers, steps, operations, elements, and/or components, but do not preclude the presence or addition of one more other features, integers, steps, operations, element components, and/or groups thereof.

The corresponding structures, materials, acts, and equivalents of all means or step plus function elements in the claims below are intended to include any structure, material, or act for performing the function in combination with other claimed elements as specifically claimed. The description of the present invention has been presented for purposes of illustration and description, but is not intended to be exhaustive or limited to the invention in the form disclosed. Many modifications and variations will be apparent to those of ordinary skill in the art without departing from the scope and spirit of the invention. The embodiment was chosen and described in order to best explain the principles of the invention and the practical application, and to enable others of ordinary skill in the art to understand the invention for various embodiments with various modifications as are suited to the particular use contemplated.

The flow diagrams depicted herein are just one example. There may be many variations to this diagram or the steps (or operations) described therein without departing from the spirit of the invention. For instance, the steps may be performed in a differing order or steps may be added, deleted or modified. All of these variations are considered a part of the claimed invention.

While the preferred embodiment to the invention had been described, it will be understood that those skilled in the art, both now and in the future, may make various improvements and enhancements which fall within the scope of the claims which follow. These claims should be construed to maintain the proper protection for the invention first described.

Claims

1. A method for providing data-driven distributionally robust optimization, the method comprising:

receiving a plurality of samples of one or more uncertain parameters for a complex system;

calculating a distribution uncertainty set for the one or more uncertain parameters;

receiving a deterministic problem model associated with the complex system that includes an objective and one or more constraints, wherein the plurality of samples is described by an unknown distribution;

creating a distributionally robust counterpart (DRC) model based on the distribution uncertainty set and the deterministic problem model;

formulating the DRC as a generalized problem of moments (GPM);

applying a semi-definite programing (SDP) relaxation to the GPM; and

generating an approximation for a globally optimal distributionally robust solution to the complex system.

2. The method of claim 1, wherein formulating the DRC as the GPM comprises:

calculating a dual minimization problem of an inner maximization problem;

transforming a feasible set of an inner minimization problem to match a structure of the feasible set of an outer minimization problem; and

reducing a minimization-minimization problem to a minimization problem, which constitutes the GPM.

3. The method of claim 1, wherein formulating the DRC as the GPM comprises:

calculating a dual-maximization problem of an inner minimization problem, transforming a feasible set of a inner maximization problem to match the structure of the feasible set of an outer maximization problem; and

reducing a maximization-maximization problem to a maximization problem, which constitutes the GPM.

4. The method of claim 1, wherein calculating the distribution uncertainty set for the one or more uncertain parameters is based on a polynomial estimate of a probability density function.

5. The method of claim 1, wherein calculating the distribution uncertainty set for the one or more uncertain parameters is based on statistical estimates for a plurality of moments of the unknown distribution of the uncertain system parameters up to an arbitrary order.

6. The method of claim 1, wherein calculating the distribution uncertainty set for the one or more uncertain parameters is based on histogram estimates for the unknown distribution of the uncertain system parameters.

7. The method of claim 1, wherein the distribution uncertainty set includes a support that is described by one or more multivariate polynomial inequality constraints.

8. The method of claim 1, wherein the objective is described as multivariate polynomial.

9. The method of claim 1, wherein the equality and/or inequality constraints are described as multivariate polynomials.

10. The method of claim 1, wherein the approximation for a distributionally robust solution includes precision level.