Applications of Stochastic Mixed-Integer Second-Order Cone Optimization

Second-order cone programming problems are a tractable subclass of convex optimization problems that can be solved using polynomial algorithms. In the last decade, stochastic second-order cone programming problems have been studied, and efficient algorithms for solving them have been developed. The mixed-integer version of these problems is a new class of interest to the optimization community and practitioners, in which certain variables are required to be integers. In this paper, we describe five applications that lead to stochastic mixedinteger second-order cone programming problems. Additionally, we present solution algorithms for solving stochastic mixed-integer second-order cone programming using cuts and relaxations by combining existing algorithms for stochastic second-order cone programming with extensions of mixed-integer second-order cone programming. The applications, which are the focus of this paper, include facility location, portfolio optimization, uncapacitated inventory, battery swapping stations, and berth allocation planning. Considering the fact that mixed-integer programs are usually known to be NP-hard, bringing applications to the surface can detect tractable special cases and inspire for further algorithmic improvements in the future.


I. INTRODUCTION
C HALLENGES, restrictions, and affects such as uncertainty [1]- [3], integrality [4]- [7], and conicity [8]- [11] arise naturally in real-world applications. For example, the (deterministic) mixed-integer secondorder cone programming (DMISOCP) models pre-sented in [12] (see also [13]- [15]) have proved to be useful in dealing with a variety of applications that involve integrality and conicity. For another example, stochastic mixed-integer linear programming (SMILP) [16] has been demonstrated to be effective in many applications involving integrality and uncertainty (see also [17]- [19] and the references contained therein). For a third example, the stochastic second-order cone programming (SSOCP) models described in [20] and the stochastic semidefinite programming models described in [21] have proved to be useful in dealing with uncertainty and conicity in many applications (see also [22] and [23]). In this paper, we show that uncertainty, integrality, and conicity can all naturally coexist together in the context of what can be termed as the stochastic mixedinteger second-order cone programming (SMISOCP).
The second-order cone is defined to contain all vectors whose first component must be at least as large as the Euclidean norm of the subvector of the remaining components. Second-order cone programming prob-lems are a tractable subclass of convex optimization problems in which the objective function is optimized over the intersection of a linear affine space and the Cartesian product of second-order cones. Interior-point algorithms are considered one of the most efficient methods developed for solving this subclass of convex optimization problems in polynomial-time (see [24]- [26], for example). The mixed-integer version of these problems, acronymized above as DMISOCP, requires some variables to be restricted to integers. Although it is NP-hard in general, DMISOCP can be quickly solved by applying algorithmic techniques that use scenario-based cuts and relaxations (see [12], for example). DMISOCP was introduced to make optimal decisions in applications with certainty in information.
In some applications, we cannot completely determine the modal since it is dependent on some information that is not available at the moment of making optimal decisions. An extra dimension of difficulty is added if we allow uncertainty in data defining DMISOCP problems. This yields the stochastic version of these problems, acronymized above as SMISOCP, which is algorithmically much less mature compared with SSOCP, SMILP, and even DMISOCP. We emphasize that there are three dimensions of difficulty in dealing with SMISOCP: Uncertainty, integrality, and conicity. We believe that these three difficulty dimensions have not been previously combined together in one class. That is why it is our belief that the SMISOCP class has not been studied yet.
In the present paper, we will demonstrate how SMISOCP applies to a variety of application areas by describing five applications leading to two-stage SMISOCP models. Specifically, we present a stochastic discrete facility location problem, a portfolio optimization with CVaR and diversification constraints, a stochastic joint uncapacitated location-inventory problem, an optimal infrastructure problem for electric vehicles with battery swap technology, and an optimal random berth allocation problem with uncertain handling time.
Since this paper mainly introduces SMISOCP from a practical point of view, we do not provide each application model with a specific algorithm to solve it. Nevertheless, a solution method for the generic SMISOCP is presented in this paper to combine existing efficient barrier [27]- [29], homogeneous [30], and infeasible [31] algorithms for SSOCP with extensions of DMISOCP algorithms [32]- [34] (see also [35]- [37]). Related aspects are tackled in an algorithmic companion paper [38]. It is our belief that bringing the applications proposed in this research to the surface can detect tractable special cases for extra algorithmic improvements in the future. This paper is structured as follows. In Section II, we introduce some notations followed by the definition of the two-stage SMISOCP problem with recourse. Section III contains the main results of the paper. Here, we describe different applications and their formulations as SMISOCP models. In Section IV, we first give a brief review of the solution algorithms available for solving SSOCP, then we present a solution method for solving SMISOCP. Section V provides some concluding remarks.

II. FORMULATION
In this part, we start by introducing certain notations, then we define the two-stage SMISOCP problem.

A. NOTATIONS AND BASICS
We utilize ";" and "," for adjoining matrices and vectors in a column and a row, respectively. As an illustration, if x and y are vectors, the following expressions are equivalent: (x T , y T ) T = (x; y).
For each vector x ∈ R n indexed from 0, we denotē x for the sub-vector containing components 1 through n − 1; therefore x = (x 0 ;x).
The nth-dimensional second-order cone (also known as the Lorentz or quadratic cone) is defined as where E n = R × R n−1 and ∥ · ∥ denotes the Euclidean norm.
The cone Q n is closed, pointed (i.e., it does not contain a pair of opposite nonzero vectors) and convex with nonempty interior in R n (see Figure 1). It is known that Q n is self-dual (i.e., it equals its dual cone). It is also known that the automorphism group of Q n acts transitively on its interior, i.e., for each u and v in the interior of Q n , there exists an invertible linear map φ : E n → E n such that φ(Q n ) = Q n (so φ is an automorphism of Q n ) and φ(u) = v. Given these properties, it is clear that Q n is a symmetric cone (see [39]). If x ∈ R n and n is known from the context, we denote x ⪰ 0 to mean that x ∈ Q n . If x ∈ R k 1 × R k 2 × · · · × R k r and k 1 , k 2 , . . . , k r are positive integers known from the context, we denote x ⪰ r 0 to mean that x ∈ Q k 1 × Q k 2 × · · · × Q k r . Therefore, x ⪰ r 0 if and only if x is partitioned conformally as x = (x 1 ; x 2 ; . . . ; x r ) and x i ⪰ 0 for each i = 1, 2, . . . , r. We also write x ⪰ r y or y ⪯ r x to mean that x − y ⪰ r 0.

VOLUME x, 2021
This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/ACCESS.2021.3139915, IEEE Access As a tool for algorithmic development, in Section IV we will use the Jordan multiplication • : E n −→ E n , which is defined as x • y = (x 0 ;x) • (y 0 ;ȳ) := x T y; x 0ȳ + y 0x . (1) In building the model of Subsection III-E, we will make use of the nth-dimensional rotated quadratic cone, which is given as [24] R n := x = (x 0 ; x 1 ;x) ∈ R × R × R n−2 : x 0 x 1 ≥ ∥x∥ 2 , A hyperbolic constraint on x ∈ R n is one that satisfies the inequality x 0 x 1 ≥ ∥x∥ 2 . The cone R n is derived by rotating the second-order cone Q n in the x 0 x 1 -plane by forty-five degrees. Note that (2) This means that the hyperbolic constraint x 0 x 1 ≥ ∥x∥ 2 is equivalent to the second-order cone constraint (x 0 + x 1 ; x 0 − x 1 ; 2x) ⪰ 0.

B. SMISOCP PROBLEM DEFINITION
Let n, m, k, l, p and q be positive integers. Let also k 1 , k 2 , . . . , k r and l 1 , l 2 , . . . , l s be positive integers such that k = r i=1 k i and l = s j=1 l j . An SMISOCP with recourse is described using deterministic data A ∈ R k×n , b ∈ R k and c ∈ R n , as well as random data T ∈ R l×n , W ∈ R l×m , h ∈ R l and d ∈ R m . The realization of this data depend on an underlying outcome ω in an event space Ω with a known probability function P. Given this data, the twostage SMISOCP with recourse is the problem where Q(x, ω) is the minimum value of the problem and Here, x and y are the decision variables for the first-and second-stage, respectively.
Due to the integrality constraints on x and y, the SMISOCP problem (3,4) is nonconvex and NP-hard in general. Below we list three special cases of this optimization problem. See also Figure 2 which visualizes the conceptual relationships between the SMISOCP problem and its special cases. • If p = q = 0, then no integrality constraints are imposed on x and y. In this case, the SMISOCP problem (3,4) is reduced to the two-stage SSOCP problem with recourse [20]. SSOCP is a convex optimization problem that can be solved in several methods in polynomial time (see [28]- [31]). • If k = r and l = s, then the conic constraints Ax ⪰ r b and W(ω)y ⪰ s h(ω) − T(w)x are reduced to linear constraints Ax ≥ b and W(ω)y ≥ h(ω) − T(w)x.
In fact, Ax ⪰ k b if and only if the vector Ax − b is nonnegative, i.e., all its components belong to the first-dimensional second-order cone Q 1 := {t ∈ R : t ≥ 0}. Similarly, W(ω)y ⪰ l h(ω) − T(w)x if and only if W(ω)y + T(w)x − h(ω) ≥ 0. Therefore, the SMISOCP problem (3,4) is reduced to the two-stage SMILP problems with recourse in this case [16], [40]- [43]. In the literature, this class of optimization problems has been studied widely [44]- [49]. Solution methods for solving SMILP can also be found in [50]- [53]. • If all data are known with certainty, the secondstage problem disappears. In this case, the SMISOCP problem (3,4) reduces to the DMISOCP problem [12]. Solution methods for solving DMISOCP can be found in [12], [32]- [34] and the references contained therein. In particular, cutting planes are one of the most successful ways to deal with the integrality constraints in this optimization problem. Solution methods available in [12], [32], [33], [35] use scenario-based cuts for DMISOCP. The results in [34] show that branch and bound is also suited for solving DMISOCP. Additionally, if the integrality constraints on x and y are restricted to be binary, we have x ∈ {0, 1} p × R n−p and y ∈ {0, 1} q × R m−q . In this case, the SMISOCP problem (3, 4) may be termed as a two-stage stochastic mixed-binary second-order cone programming (SMBSOCP) problem with recourse. SMBSOCP is known to be as "hard" as the SMISOCP problem.
We refer the reader to two papers that deal with important special cases of or related to the SMISOCP problem (3,4): The first paper is by Luo and Mehrotra [54] which extends the work of Sen and Sherali [51] for SMILP and proposes a decomposition method for (3,4) in which x ∈ {0, 1} p × R n−p in the first-stage problem in (3,4) and y ∈ Z q × R m−q in the second-stage problem in (3,4). The second paper is by Bansal and Zhang [55] which proposes scenario-based cuts for (3,4) in which n = p and k = r in the first-stage problem in (3,4) and the l 2 -norms are generalized to l p -norms in the second-stage problem in (3,4).

III. APPLICATIONS
Each of the following subsections will describe an application accompanied with SMISOCP model formulation.

A. STOCHASTIC DISCRETE FACILITY LOCATION PROBLEM
In a facility location problem (FLP), we are interested in choosing a location or locations to construct a new facility or new facilities so that desired distances from the new facility(ies) to existing facilities are minimized. In a Euclidean FLP, the desired distances are Euclidean distances. FLP arises in many applications including, but not limited to, locating regional campuses, educational sites, public health clinics, pharmacies and hospitals, placing servers in probabilistic networks and data centers, establishing and co-locating wireless telecommunication facilities, clustering data of high dimensionality, performing any type of location analysis for placement of new facilities of varying sizes and densities, and so on.
FLP arises also in many applications where the uncertainty is an important challenge. To realize this more, we describe a concrete example of the stochastic version of this application. We consider how to distribute new appropriate facilities among COVID-19 healthcare existing facilities across a small country, like Jordan, keeping in mind that COVID-19 is very unpredictable and easily spread. Given the pandemic circumstances, the healthcare existing facilities are either fixed or random. The fixed facilities are the constructed clinics or hospitals in the Jordanian cities. The random facilities are movable healthcare facilities equipped with the necessary medical supplies which enable patients to be treated immediately. As shown in Figure 3, the movable healthcare facilities are random existing facilities that can be moved from their sites to other sites depending on the prediction of the COVID-19 spreading rate in a model country. In the stochastic COVID-19 healthcare facility location problem, we know the specific locations of the fixed healthcare facilities, but we do not know the specific locations of the random healthcare facilities at the time of formulation. Instead, we know the realizations of the random facilities with a predefined probability p ∈ (0, 1).
The location of a new facility can be a specific location for dropping off supplies, staff housing units, a healthcare unit with experts in treating mood and anxiety disorders, a fully functional advanced comprehensive unit, etc. In order to guarantee accuracy, the practitioners will use only the most recent information on the community's COVID-19 epidemic outbreak. This may require increasing the initial number of random facilities. This increase, based on the number of fixed facilities and the initial number of random facilities, results in an additional cost.
Stochastic FLP arises in many applications where the integrality is an important restriction. Below we describe another, but discrete, concrete example of this generic application. In a discrete facility location problem with a risk of disruption, we specify a set of demand points (customers) to service and a set of potential facility points (or sites) to serve them. At the time of formulation, not all the potential facility points have constructed facilities, nor have opened facilities. Indeed, some constructed facilities are unavailable or closed due to disruptions caused by some factors such as weather, labor strikes, pandemic outbreak, natural disasters, terrorist attacks, etc. The discretion lies in the determination of whether a potential facility point has a constructed facility. The uncertainty lies in the determination of whether some of the constructed facilities are open. We are required to pick a subset of potential sites for facility construction and to choose a subset of facilities to open so that the sum of the distances between any demand point and its nearest opened facility, and the sum of the opening costs of the selected facilities are both minimized. Figure 4 shows a 49-demand-points solution to this problem in an interactive USA map. In this figure, the constructed facilities that are open with uncertainty have randomness in terms of openness. At the present time, the realizations of such facilities are unknown, but they will become known at some point in the future. The model contains these realizations with a prescribed probability p ∈ (0, 1).
In addition to the above two concrete examples, many other examples can also be referenced. Generally speaking, the locations of some existing facilities are not totally specified in a stochastic FLP. Such unspecified locations are dependent upon information that was not available at the time of making the decision and will not become available until a later time. In a discrete FLP, some or all new facilities must not be placed in any location, but must be located in specific locations, called nodes, in the solution space. The (deterministic) Euclidean FLP is often presented as an application of (deterministic) second-order cone programming (see [56], for example). The stochastic Euclidean FLP was presented in [20] as an application of SSOCP. In this subsection, we present the stochastic discrete version of the generic Euclidean FLP as an application of SMISOCP.
In summary, let a 1 , a 2 , . . . , a r be fixed points in R n symbolizing the coordinates of r existing fixed facilities, andã 1 (ω),ã 2 (ω), . . . ,ã s (ω) be random points in R n that represent the coordinates of s random facilities on which their realizations are dependent on an underlying outcome ω in an event space Ω with a known probability function P. Let us assume that right now the realizations of s random facilities are not known, and over the course of time their realizations will be known.
The goal is to add m new facilities, x 1 , x 2 , . . . , x m ∈ R n , that minimize each of the following sums: • The weighted sums of the Euclidean distance between each one of the new facilities, x j for j = 1, 2, . . . , m, and each one of the existing fixed facilities, a i for i = 1, 2, . . . , r. This requirement can be met by adding the inequalities t i j ≥ ∥x j − a i ∥ for i = 1, 2, . . . , r and j = 1, 2, . . . , m. Equivalently, using our notations, this is satisfied by adding the constraints t 1j ; x j − a 1 ; . . . ; t r j ; x j − a r ⪰ r 0, for j = 1, 2, . . . , m, where the double summation m j=1 r i=1 w i j t i j is minimized. The parameter w i j specifies the weight of the ith existing facility and the jth new facility, for i = 1, 2, . . . , r and j = 1, 2, . . . , m. We point out that if no interaction exists among the new facilities, we omit these constraints from our model because this requirement is no longer present.
Notably, our decision should be made prior to the availability of realizations of the s random facilities. Only the most recent information concerning random facilities is utilized to ensure accuracy. This may need an increasing or decreasing number of new facilities depending on the availability of the most recent data on random facilities. For simplicity, let us assume that the number of new facilities was previously determined and fixed (which is m).
One variant of this application is to require, in addition to the above three requirements, each new facility x j to be placed at specific locations within the solution space. Such specific locations are represented by fixed points v for each j = 1, 2, . . . , m. This leads to VOLUME x, 2021 the following SMBSOCP model.  (1) ; z (2) ; . . . ; z (m) ; x 1 ;  (1) ; y (2) ; . . . ; y (m) ; x 1 ; Another variant of this application is to require some of the new facilities to attain integer values, i.e., x j has to be integer-valued for each j ∈ ∆ ⊂ J. Without loss of generality, let us assume that x 1 , x 2 , . . . , x |∆| ∈ Z n . In this case, we are interested in the following model.
As mentioned in Section II, a solution for a model such as that in (5) and (6) is now possible. Specifically, Model (5,6) can be solved by the algorithm proposed in [54]. The SMISOCP model (7,8) is also intractable for an extensive form formulation (see the form (36)).

B. PORTFOLIO OPTIMIZATION WITH CVAR AND DIVERSIFICATION CONSTRAINTS
Portfolio optimization is one of the most important applications of second-order cone programming. One of its paradigms is maximization of the expected return over a single period subject to conditional value-at-risk (CVaR) constraints. These constraints are introduced to limit the downside risk of the tracking portfolio. In our application model, we also consider diversification-bysector constraints, which are used to limit the minimum amount invested in each asset from each economic sector, and buy-in threshold constraints, which are used to prevent the investor from investing very small amounts in a single asset. The last two types of constraints mentioned above were considered in [57].
Portfolio optimization is frequently framed as an application of (deterministic) second-order cone programming (see [35], for example). The stochastic perspective of this application was presented in [20] as an application of SSOCP by maximizing expected returns over two periods or stages. The discrete version of this application was presented in [12] as an application of DMISOCP by adding diversification-bysector constraints. The application model presented in this subsection combines both perspectives by adding diversification-by-sector constraints and buy-in threshold constraints over two stages, which results as an SMISOCP model.
We consider cash (index 0) and n risky assets from L distinct sectors in our portfolio during a two-stage period. Let x + , x − ∈ R n+1 be the buy and sell transactions, respectively, throughout the first stage. We utilize the following parameters over the first stage: we use ξ ∈ {0, 1} L to indicate if each sector has sufficient investments, use δ ∈ {0, 1} n for buy-in threshold constraints, use u ∈ R n+1 to denote the current portfolio holdings, and let p ∈ R n+1 denote the expected rates of return. For simplicity, we also assume that p is Gaussian with known meanp and covariance Σ, so the return over this stage is the Gaussian random variable r = p T (u + ∆x) with meanr =p T (u + ∆x) and variance σ = (u+∆x) T Σ(u+∆x), where the difference ∆x := x + −x − is the trading shares generate cash in the first stage.
Let also y + , y − ∈ R n+1 be the buy and sell transactions, respectively, throughout the second stage. We use the following parameters over the second stage: We use ζ ∈ {0, 1} L to denote whether if each sector has sufficient investments, use η ∈ {0, 1} n for buyin threshold constraints, use v ∈ R n+1 to denote the current portfolio holdings, and let q(ω) ∈ R n+1 denote the random expected rates of return whose realization depends on an underlying outcome ω in an event space Ω with known probability function P. For simplicity, we also let q(ω) be Gaussian with known meanq(ω) and covarianceΣ(ω), so the return over this stage is the Gaussian random variabler(ω) = q T (ω)(v + ∆y)

The first-stage
The second-stage · Expected rates of return p.
Portfolios will be kept until t 3 after they were determined at t 1 and potentially revised at time t 2 . t 2 t 3 time Figure 5: Two-stage portfolio optimization with CVaR constraints.
with meanr(ω) =q T (ω)(v + ∆y) and varianceσ(ω) = (v+∆y) TΣ (ω)(v+∆y), where the difference ∆y := y + −y − is the trading shares generate cash in the second stage. We penalize the buy and sell transactions throughout the two stages by quadratic convex transaction costs of Garleanu and Pedersen [58]. The transaction costs associated with the first and second stages are given by ∆x T Λ ∆x and ∆y TΛ (ω) ∆y, respectively, where Λ,Λ(ω) ∈ R (n+1)×(n+1) are symmetric positive definite matrices measuring the level of trading costs. The matrix Λ is obtained as a positive multiple of Σ, and Λ(ω) is obtained as a positive multiple ofΣ(ω).
Assume that we do not know the Gaussian vector's realization q(ω) at the present time, but that it will be known at a certain point in time in the future.
Our objective is to select the optimal trading strategy to maximize the expected total return at the end-ofsecond-stage. Our choices of portfolios involve the Markowitz trade-off between their expected return and investment market risk (or volatility) [59].
From the above discussion, it can be seen that the objective function of the first-stage problem is where Q(∆x, ω) is the maximum value of the secondstage problem, which has the objective function It is beneficial to provide continuous relaxations of the above objective functions that replace their quadratic forms by linear forms. This can be achieved via including only linear and second order cone constraints as is shown below: Similarly, the objective function of the second-stage problem can be stated as To avoid falling below the expected wealth levels, we need to add M CVaR constraints at the end of the first stage, and N CVaR constraints at the end of the second stage. For i = 1, 2, . . . , M and j = 1, 2, . . . , N, let α i and α j represent unwanted expected wealth levels at the end of the first and second stages, respectively, and let β i andβ j denote the first and second stage's minimal probability, respectively. Assuming the above data is given, our goal is to maximize the end-of-second-stage expected total return less transaction costs. In order to limit our risk, each CVaR constraint i, i = 1, 2, . . . , M, asks that our expected wealth at the end of the first stage is above α i with a probability of at least β i , and each CVaR constraint j, j = 1, 2, . . . , N, asks that our expected wealth at the end of the second stage is abovẽ α j with a probability of at leastβ j .
Given the above data, we need to determine the decision variables x + and x − so that the CVaR constraints are satisfied at the end of the first stage, and determine the decision variables y + and y − so that the CVaR constraints are satisfied at the end of the second stage. These determinations need to be made before the realizations of the random expected rates of return become available.
To be more precise, assume that 0 < t 1 < t 2 < t 3 are given times. Let t 1 and t 2 be the starting times of the first and second stages, respectively, and t 2 and t 3 be the ending times of these two stages, respectively, as shown in Figure 5. At time t 1 , the investor is given the vector of returns p and is required to determine the portfolio vectors x ± . At time t 2 , the investor can revise some, if not all, entries of the portfolio vectors x ± to get the portfolio vectors y ± . The revised portfolios y ± are kept until time t 3 . The objective is the determination of the portfolios at VOLUME x, 2021 time t 1 in anticipation of the revision of some or all of them at time t 2 . The CVaR constraints (9) are equivalent to the secondorder cone constraint is the cumulative distribution function for a standard normal random variable. A common proof of the above reformulation of the CVaR constraints (9) is as follows (see also [12], [60]): We can write the CVaR constraint P(r ≥ α i ) ≥ β i as Since (r−r)/ √ σ is a zero mean unit variance Gaussian random variable, the above probability is simply Using the symmetry of the standard distribution function and in light of the equality √ σ = ∥Σ 1/2 (u + ∆x)∥, the conicality constraint (11) follows from the following equivalences, for i = 1, 2, . . . , M.
Similarly, one can show that the CVaR constraints (10) are equivalent to the conicality constraint It is now time to introduce the previously mentioned diversification constraints, which force the investor to diversify its portfolios at each stage by buying assets in at least L min different sectors. In the first stage, we link every asset i with a sector k such that the sets S k , 1 ≤ k ≤ L, of assets affiliated with a sector k produce an exact partition of {1, 2, . . . , n}. We also associate a variable ξ k ∈ {0, 1} with each sector k, so it is equal to 1 if and only if s min is a minimum predefined level of the amount of sector k over the first stage. That is, where 1 S k ∈ {0, 1} n+1 is the indicator vector of S k whose jth entry has value 1 if j ∈ S k and 0 otherwise. Additionally, our portfolio's total allocation of assets to an economic sector k cannot be greater than one. As a constraint involving a binary variable ξ k , this can be formulated as In the second stage, any asset j is associated with a sector k as for those of the first stage, and a variable ζ k ∈ {0, 1} is associated with each sector k as follows: wheres min (ω) is a minimum predefined level of the amount of sector k over the second stage, whose realization depends on an underlying outcome ω in an event space Ω with known probability function P.
Similarly, as a constraint involving a binary variable ζ k , this can be formulated as In order to satisfy the diversification requirements to detain "representative" positions in at least L min sectors, we should additionally add the cardinality constraints where 1 is the vector of ones of an appropriate dimension.
We now consider the buy-in threshold constraints. As mentioned at the beginning of this subsection, these constraints are introduced to avoid the investor from investing very small amounts in a single asset at each stage. As constraints involving binary variables δ and η, these requirements can be formulated as where u min and v min are predetermined proportion of the available capital budgets to be invested.
Our portfolio model is based on the simplest assumption that we require full portfolio allocations in the stocks at the end of each stage, i.e., 1 T (u + ∆x) = n j=0 (u j + ∆x j ) = 1 and 1 T (v + ∆y) = n j=0 (v j + ∆y j ) = 1. Additionally, for our portfolio vectors, we also assume that we can take limited short positions for each nonliquid asset to allow for short sales of them and enhance returns during volatility. This adds the constraints u + ∆x ≥ −s and v + ∆y ≥ −s, where s ands represent the short position limits at the first and second stages, respectively, for each nonliquid asset. Finally, we need the variables associated with the buy and sell transactions to be nonnegative at each stage, i.e., x ± ≥ 0 and y ± ≥ 0.
Given the above data and requirements, the stochastic portfolio optimization problem with CVaR and diversification constraints over two stages can be cast as the  Figure 6: A two-stage stochastic inventory distribution process.

Stage II
following SMBSOCP model.
As indicated earlier, since Model (12,13) has binary variables, a solution for such a model is now possible by applying the decomposition algorithm developed in [54].

C. STOCHASTIC JOINT UNCAPACITATED LOCATION-INVENTORY PROBLEM
Inventory management is important in numerous industries such as E-commerce, consumer goods, car renting, food distribution, or any industry or trade that deals with inventories. We consider the basic joint uncapacitated location-inventory problem over two periods (or stages) with some uncertainty in data relating to retailers. More specifically, we consider a supply chain in which assets are shipped by suppliers from distribution centers to fixed retail stores over Stage I, and are also shipped to fixed and random retail stores over Stage II (see Figure 6). A given period or stage can be a few months or a year for example. Shen et al. [61] and Daskin et al. [62] studied the same problem over one period with nonlinearity arising from safety stock costs. Atamturk et. al [63] were the first to formulate the one-stage joint location-inventory model proposed in [62] and [61] as a DMISOCP model.
Our problem is the following. We are given a number of candidate locations for opening distribution centers. We are also given a number of existing retailers over Stage I and a number of existing and/or random retailers over Stage II, each with uncertain product demand. At each stage, it is necessary to decide the number of distribution centers that need to be opened, the kind of retailers to be allocated to each center, in addition to the required level of safety stock to be maintained in order to minimize total inventory costs, opening location, and shipment while providing a certain level of service. In each stage in the model, we assume that the shipments are made directly from distribution centers to retailers, and that the demand is independent and Gaussian at each retailer.
In building the model, the distributional constraint that we consider in each stage is that each retail store is supplied from exactly one open distribution center. To model this constraint, let I be the set of fixed existing retailers in the first stage, andĨ ⊆ I be the set of fixed existing retailers in the second stage. We also let J be the VOLUME x, 2021 if location j is opened for a new distribution center, 0, otherwise; if a fixed retailer i is assigned to a new distribution center located at location j, 0, otherwise; if a random retailer k is assigned to a new distribution center located at location j, 0, otherwise. (14) set of candidate locations for opening distribution centers, and K be the set of random retailers in the second stage whose realization depends on an underlying outcome ω in an event space Ω with known probability function P. We introduce the decision variables x, y (1) , . . . , y (|I|) , z (1) , . . . , z (|K|) ∈ {0, 1} |J| , which are defined in (14). Given the above decision variables, the requirement that each retailer is supplied from exactly one distribution center and the requirement that retailers are only assigned to open distribution centers are modeled in Stage I as j∈J y (i) j = 1, and y (i) j ≤ x j , i ∈ I, j ∈ J, respectively, and are modeled in Stage II as j∈J y In Stage I, we minimize the following costs: • The fixed cost of locating the jth distribution center for j ∈ J. This objective term can be written as f j x j , where f j is the monthly fixed cost of locating a new distribution center at location j. • The cost of shipping from the central plant to the jth distribution center for j ∈ J. a j represents the unit cost of shipment from the central plant to the jth distribution center, µ i is the mean of daily demand at the ith retailer, and β is the weight associated with the inventory costs. After that, this objective term is written as β i∈I a j µ i y (i) j .
• The cost of shipping from the jth distribution center to retailers, for j ∈ J. Let d (i) j be the unit cost of shipment from the jth distribution center to the ith retailer, then this objective term is written as • The expected working inventory cost at the jth distribution center for j ∈ J. Following Shen et al. [61], we describe the process by which this objective term is derived. The expected working inventory cost includes the total shipment cost per year, the fixed cost of placing n orders per month, and the average working inventory cost. This objective term is written as j is the monthly expected demand at the jth distribution center, F j is the fixed cost of placing an order at the jth distribution center per month, v j (x) is the shipment cost function for x units from the central plant to the jth distribution center, h is the unit inventory holding cost per month, and θ is the weight associated with the transportation costs. Assume that v j (·) is linear in its parameters and is specified as v j (x) = a j x + g j , where g j is the fixed cost per shipment from the central plant to the jth distribution center. Then we have the expression: Taking the derivative of the above expression with respect to n, we obtain F j + βg j − θhD j /2n 2 . Equating this to zero and solving for n, we get n = θhD j /2(F j + βg j ). By substituting n into (15), we get 2θhD j (F j + βg j ) + βa j D j . Note that βa j D j = β i∈I a j µ i y (i) j already exists in the objective. Therefore, this objective term is reduced to • The expected safety stock cost at the jth distribution center for j ∈ J. Let σ i be the standard deviation of daily demand at the ith retailer, z α be the standard normal deviation associated with service level α, and L j be the lead time in days at the jth distribution center. Then this objective term is written as Similarly, in Stage II, the following parameters are used for i ∈Ĩ, k ∈ K and j ∈ J: f j (ω): the monthly fixed cost of locating a new distribution center at location j, a j (ω): the unit cost of shipment from the central plant to the jth distribution center, µ i : the mean of daily demand at the ith fixed retailer, µ k (ω): the expected mean of daily demand at the kth random retailer, β(ω): the weight associated with the inventory costs, d (i) j : the expected unit cost of shipment from the jth distribution center to the ith fixed retailer, d (k) j (ω): the unit cost of shipment from the jth distribution center to the kth random retailer, F j (ω): the fixed cost of placing an order at the jth distribution center per month, h(ω): the unit inventory holding cost per month, θ(ω): the weight associated with the transportation costs, g j (ω): the fixed cost per shipment from the central plant to the jth distribution center, σ i : the standard deviation of daily demand at the ith fixed retailer, σ k (ω): the expected standard deviation of daily demand at the kth random retailer, z α (ω): the standard normal deviation associated with service level α, L j (ω): the lead time in days at the jth distribution center, ; the expected monthly demand at the jth distribution center, v j (x, ω) :=ã j (ω)x +g j (ω); a shipment cost function for x units from the central plant to distribution center j.
The realizations of the above unspecified parameters depend on an underlying outcome ω in an event space Ω with a known probability function P. Such undetermined parameters are dependent on information that was not available at the time of making the decision, and will become available at a later point in time. At the beginning of Stage II, the investor can revise some, if not all of these parameters. The revised parameters are kept until the end of Stage II. Likewise, in Stage II the following costs, which can be derived as those costs in Stage I, are minimized for i ∈Ĩ, k ∈ K and j ∈ J: • The expected fixed cost of locating the jth distribution center:f j (ω)x j . • The expected cost of shipping from the central plant to the jth distribution center: • The expected cost of shipping from the jth distribution center to (fixed and random) retailers: • The expected working inventory cost at the jth distribution center: j .
• The expected safety stock cost at the jth distribution center:zα The objective is to minimize these costs at the beginning of Stage I in anticipation of the revision of some or all of the unspecified parameters at the beginning of Stage II. Let alsô The description of the two-stage stochastic joint location-inventory problem and the mathematical modeling of its constraints and objective shown above lead to the following two-stage stochastic nonlinear program.
The nonlinear terms in (16) and (17) can be handled VOLUME x, 2021 by introducing the auxiliary variables s j ,s j , t j ,t j ≥ 0 and using the fact that (y (i) j for i ∈ I and k ∈ K. Given this, Problem (16,17) can be cast as the following two-stage SMBSOCP model.
where Q(x; y (1) ; . . . ; y (|Ĩ|) , ω) is the minimum value of the problem We emphasize that Model (18,19) accepts only binary variables and is intractable for an extensive form formulation. Therefore, a solution for such a model is now possible by applying the decomposition algorithm developed in [54].

D. OPTIMAL INFRASTRUCTURE PROBLEM FOR ELECTRIC VEHICLES WITH BATTERY SWAP TECHNOLOGY
Electric vehicles are today an attractive option for many car shoppers and very affordable to them, but the operation of charging the electric vehicles is still the most significant challenge faced by drivers because it requires a significant amount of time. To overcome this obstacle, a new infrastructure strategy has been developed recently, seeking to exchange depleted batteries in electric vehicles with full ones in the middle of long trips. In this part, we describe the problem of creating infrastructure, offering service for electric vehicles with battery swap mode, and providing coverage for stations with certain and uncertain arrival rates over two stages.
We point out that the corresponding single-stage problem was described in [64] as an application of DMISOCP. The objective is to build a number of battery swapping stations at strategic locations along an existing freeway network.
We consider an existing freeway network connecting a number of cities. Let P be the set of intercity travel paths on the network, and Q be the set of subpaths of paths p ∈ P and length longer than d, where d is the maximum travel distance allowed with one full vehicle charge. Any electric vehicle traveling a long distance needs to get access to a swap station before its battery runs out. Therefore, it must visit at least one station along each subpath of a travel path that is longer than a distance of d/2 in a round trip, or longer than a distance of d in a one-way trip. Figure 7 shows a highway between two Jordanian cities that consists of 10 linked segments, which are defined between 9 adjacent exits, and 6 corresponding subpaths of the path. The first subpath consisting of segments 1-5, and the next subpaths consisting of segments 2-6, 3-7, 4-8, 5-9, and 6-10. Each one of these subpaths must pass through at least one swapping station.
The number of electric vehicles traveling along the network is known along some portion of paths and is uncertain along other portions of paths. Let J be the set of candidate locations for swapping stations with known and fixed demand, andJ be the set of candidate locations for swapping stations with uncertain demand. We consider a problem over two stages as follows. In Stage I, the values of the arrival rates at swapping stations at all sites j ∈ J are precisely known, but such values at swapping stations at all sites j ∈J are not precisely known. Therefore, all swapping stations at sites j ∈ J are built in Stage I. At the beginning of Stage II, the realizations of the uncertain arrival rates are precisely observed, hence we can revise some information about swapping stations under uncertainty to build all stations and meet demand. Figure 8 shows a concrete example of this application. We introduce the decision variables x, y (p) , z (q) ∈ {0, 1} |J|+|J| , for p ∈ P and q ∈ Q, which are defined in (20) for j ∈ J ∪J.
We assume that the swap mode of electric vehicles traveling along the path p ∈ P enter the freeway network via a Poisson process with a rate λ (p) . Accordingly, at a swapping station at site j ∈J, electric vehicles arrive and request for swapping service according to a Poisson process with rate λ j = p∈P λ (p) y (p) j which is a precisely known quantity now in Stage II. The arrival rates λ (p) , p ∈ P, depend on a number of mutually independent random factors, ξ l (ω), l = 1, . . . , L, known as the primitive uncertainties. In Stage I, only the means (denoted by µ l , l = 1, . . . , L), covariance matrix (denoted by Σ with the individual standard deviations σ l , l = 1, . . . , L), and supports of the element of the random vector ξ(ω) = (ξ 1 (ω), . . . , ξ L (ω)) are known, but not the precise distribution. At the beginning of Stage II, the realization of ξ(ω) is precisely observed.
We also use the following parameters for j ∈ J ∪J, p ∈ P and q ∈ Q: f j : the annualized fixed cost incurred if location j is opened for a swapping station, h: the annualized holding cost per battery, a (q) j : a binary parameter indicating whether a station at site j is along the qth subpath (= 1) or not (= 0), b p,q : a binary parameter indicating whether the qth subpath is part of the pth path (= 1) or not (= 0), I j (y): the number of batteries to be stocked at a station at site j given the assignments of stations to paths, g j : the maximum number of batteries that can be safely accommodated by the electric grid at site j.
In stage I, we assume that the quantities I j (y) are known and fixed for j ∈ J, but they are random for j ∈J and their realizations depend on an underlying outcome ω in an event space Ω with a known probability function P. In other words, the quantities I j depend on both the choice locations and demand realizations when j ∈ J, and do not depend on the chosen of locations nor on demand realizations when j ∈J. This is represented as follows: I j (y) = I j for j ∈ J, and I j (y) = I j (y, ω) for j ∈J. Specifically, in Stage I, the charging service provider knows the precise values of the arrival rates at swapping stations at all sites j ∈ J, but does not know such values precisely at swapping stations at all sites j ∈J. After constructing all the swapping stations at sites j ∈ J in Stage I and precisely observing the realizations of the uncertainties at the beginning of Stage II, the service provider can revise some, if not all of the quantities I j (y, ω), j ∈J, in order to stock batteries at swapping stations at all sites and meet demand.
In building the model, we make sure that the following constraints are satisfied: • Electric vehicles traveling any subpath q ∈ Q need to visit at least one swapping station. This can be modeled as j∈J∪J a (q) j z (q) j ≥ 1 for each q ∈ Q. • If a subpath q ∈ Q, with a swapping station at site j, is a part of multiple paths p ∈ P, then each of those paths inherit the station at j along q. This can be modeled as y j for j ∈ J ∪J, q ∈ Q, p ∈ P. • Electric vehicles are assigned to open stations only.
This can be modeled as y (p) j ≤ x j for j ∈ J ∪J, p ∈ P. • The number of batteries at a station j ∈ J is bounded by the maximum allowable number of simultaneous parallel recharges permitted by the electric grid. This can be modeled as I j ≤ g j for j ∈ J. Dep. Des.
Battery swapping station with certain demand.
Battery swapping station with uncertain demand.

Des. Destination.
Reachable swapping stations. x j = 1, if a swapping station is located at site j, 0, otherwise; if an electric vehicle traveling along a path p will visit swapping station at site j, 0, otherwise; if an electric vehicle traveling along a subpath q will visit swapping station at site j, 0, otherwise. (20) • With a high probability of at least 1 − ϵ g , the number of batteries at a station j ∈J is bounded by the maximum allowable number of simultaneous parallel recharges permitted by the electric grid. This constraint acts on the worst-case (infimum) probability for all possible distributions in F, where F is the family of all possible distributions with the specified means and variances. This can be modeled as We consider the family F to be nonempty, i.e., the given means and variances are calculated using a univariate distribution.
The objective of the model is to minimize the following costs: • The cost of opening a swapping station at site j ∈ J ∪J and equipping it with swapping machinery. This objective term is written as f j x j for j ∈ J ∪J. • The annualized holding cost at a swapping station at site j ∈ J. This objective term is in the first-stage problem and is written as h I j for j ∈ J. • The expected annualized holding cost at a swap-ping station at site j ∈J under the worst possible (i.e., supremum) distribution over F. This objective term is in the second-stage problem and is written as h sup P∈F E P [I j (y, ω)], j ∈J.
The description of the two-stage stochastic optimal infrastructure planning problem for electric vehicles with battery swap mode and the mathematical modeling of its constraints and objective shown above lead to the following two-stage stochastic nonlinear program.
where E[Q(x; y; z, ω)] = Ω Q(x; y; z, ω)P(dω), and Q(x; y; z, ω) is the minimum value of the problem j , j ∈J, p ∈ P, q ∈ Q, y (p) j ≤ x j , j ∈J, p ∈ P, inf P∈F P P I j (y, ω) ≤ g j ≥ 1 − ϵ g , j ∈J, Our aim now is to formulate both the nonlinear term in problem (24) in regards to its objective function and the chance constraint (the fourth constraint) in (24) as second-order cone constraints. We make the following assumptions, which are quite standard [64]. Assumption 1: The path-based demand arrival rate, λ (p) , is a linear function of a number of mutually independent random factors, ξ l (ω), l = 1, 2, . . . , L. That is, l , l = 0, 1, . . . , L, are known constants. Assumption 2: In at least α (> 0.5) proportion of the swap requests, the electric vehicle picks up a battery that has been recharged for at least t time units. Assumption 1 is for convenience of obtaining tractable constraints. The requirement in Assumption 2 is equivalent to requiring that the number of Poisson demand arrivals in a period of t time units is smaller than the quantity I j (y, ω) with probability α.
We approximate the Poisson demand with a Normal distribution that matches the mean and variance. This approximation is quite standard (see for example [64]), and is accurate especially when the demand rate is large, which is the case in practice on freeway networks. To satisfy these service requirements, the number of batteries needed is approximately where Φ(·) is the cumulative distribution function of standard Normal. Therefore, for j ∈J, we have where the second and third equalities follow from the definitions of λ j and λ (p) , respectively, the fourth equality follows from the fact that y j , and the value of Ψ was obtained from Proposition 2 in [64] and is given in (22).
It immediately follows that minimizing the quantity sup P∈F E P [I j (y, ω)] is equivalent to minimizing a vari-VOLUME x, 2021 able, sayṽ j , subject to the constraint where j ∈J. Now, for j ∈J, let y j be the vector in R |P| whose pth component is y (p) j for p ∈ P, and D be the diagonal matrix in R |P|×|P| whose (pth, pth) entry is ( L l=1λ (p) l µ l ) 1/2 for p ∈ P. Then the constraint in (26) can be written as the second-order cone constraint Now, we formulate the chance constraint in (24) as a second-order cone constraint. Note that, for j ∈J, we have inf P∈F P P (I j (y, ω) ≤ g j ) = inf P∈F P P ( p∈P L l=1λ (23), whereĝ j is a constant that is given bŷ and the second and the fourth equalities in (23) follow from observing that the number of batteries needed, as specified in (25), is a strictly increasing function in the demand rate λ j . It immediately follows that the chance constraint P( p∈P L l=1λ j ) is a zero mean unit variance Gaussian random variable, the above probability is simply Φ((ĝ j − p∈P L l=1λ Because √ σ = ∥Σ 1/2 y j ∥, the given chance constraint can be written aŝ or equivalently In the second stage problem, as we used the parametersṽ j , j ∈J, we also use v j to denote the previously known quantities I j , j ∈ J, in the first-stage problem. Given the above modeling framework, Problem (21,24) can be cast as the following two-stage SMBSOCP model.
min j∈J ( f j x j + hv j ) + E Q(x; y; z, ω) where Q(x; y; z, ω) is the minimum value of the problem min j∈J ( f j x j + hṽ j ) (28) Because the constraints in the first-stage problem (27) are linear with integer variables and without the use of second-order cone relaxations, a solution for Model (27,28) is now possible. One way to approach this model is by applying scenario-based cuts proposed in [55].

E. OPTIMAL RANDOM BERTH ALLOCATION PROBLEM WITH UNCERTAIN HANDLING TIME
In recent years, container shipping has seen rapid growth especially during the pandemic of COVID-19.
For instance, Figure 9 shows that, contrary to expectations, container shipping in Shanghai has quickly bounced back from its slowdown at the beginning of the pandemic. Vessel emissions and fuel cost are today real challenges experienced in the berth allocation, and their reduction is a common target and a key priority in the port-shipping coordination. In this part, we consider this berthing allocation problem.
Du et al. [65] studied the determination of the optimal berthing positions and order for vessels waiting at a container terminal to minimize the vessels' waiting and fuel consumption. They regard the vessels handling times as previously-known (deterministic) times when formulating this problem. Due to some uncertainties, such as mechanical problems, system failures, bad weather, and route changes, the scheduled vessels  handling times may diverge from the deterministic ones, which leads to significant changes of other terminal operations and makes the baseline berthing plan no longer feasible. In this part, we consider the optimal random berth allocation problem, in which the vessels handling times and the required departure times of the vessels have known deterministic and random frames.
In the berth plan, we assume that we are given f vessels, say V 1 , . . . , V f , whose handling times, say h 1 , . . . , h f , are precisely known at the time of formulation. We also assume that we are given r vessels, say V f +1 , . . . , V f +r , whose handling times, say h f +1 (ω), . . . , h f +r (ω), are not precisely known at the time of formulation, and their realizations depends on an underlying outcome ω in an event space Ω with a known probability function P.
We presume that at present we do not know the realizations of the handling times and the realizations of the requested departure times of the vessels V f +1 , . . . , V f +r at the present, and that at some point in time in the future the realizations of these times become known. We also assume that we need to determine the leftmost berthing position of the vessels V 1 , . . . , V f , V f +1 , . . . , V f +r , and their start time of berthing given the handling times h 1 , . . . , h f and the realizations of the handling times h f +1 (ω), . . . , h f +r (ω). However, this decision must be made before the realizations of the handling times h j (ω), j = f + 1, . . . , f + r, and the realizations of the requested departure times become available. Therefore, when these realizations become available, the start times of berthing that have already been determined may or may not be optimal. At that stage, we assume that it is permitted to modify the start time of berthing of the vessels V f +1 , . . . , V f +r (but not their leftmost berthing positions), as necessary, to ensure that all berthing start times are optimal.
Define the following decision variables for i, k = 1, . . . , f, f + 1, . . . , f + r (i k): x i : the leftmost berthing position of the vessel V i , y i : the start time of berthing of the vessel V i , a i : the arrival time of the vessel V i , α ik : a binary variable indicating whether V i is positioned left of V k along the wharf (= 1) or not (= 0), β ik : a binary variable indicating whether V i is positioned below V k along the timeline (= 1) or not (= 0).
We also use the following parameters: In building the model, we make sure that the follow-ing constraints are satisfied. See also Figure 10 which visually shows an example clarifying some of these constraints graphically.
• All vessels should be berthed within the wharf's boundary. This can be written as x i + l i ≤ L for i = 1, 2, . . . , f + r. • For i, k = 1, 2, . . . , f + r, i k, if vessel V i is positioned along the wharf to the left of vessel V k , the rightmost berthing position of vessel V i must be to the left of the leftmost berthing position of vessel V k . Otherwise, the rightmost berthing position of vessel V i must be limited by the leftmost berthing position of vessel V k plus the wharf length of the container terminal. This can be modeled as for i, k = 1, 2, . . . , f + r, i k. For example, in Figure  10, it can be seen that the inequality in (29) is satisfied for the vessels V i , V k and V m as well as for the vessels V i , V j and V m , but it is not satisfied for the vessels V k and V j . • A vessel should not be berthed before its arrival time. This can be written as a i ≤ y i , for i = 1, 2, . . . , f + r. See Figure 10. As a matter of fact, the deterministic and stochastic times in Figure 10 are existing on two different axes for more focused picture, but indeed they both occur on the same timeline. • For i = 1, 2, . . . , f and k = 1, 2, . . . , f + r, i k, if vessel V i is positioned below vessel V k along the timeline, the departure time of vessel V i must precede the berthing time of vessel V k . Otherwise, the berthing time of vessel V k should not be limited by the departure time of vessel V i . Let M be an arbitrary large constant, then this can be modeled as Figure 10, it can be seen that vessels V i and V k satisfy the inequality in (30).
. . , f +r and k = 1, 2, . . . , f +r, i k, if vessel V j is positioned below vessel V k along the timeline, the expected departure time of vessel V j must precede the berthing time of vessel V k . Otherwise, the expected departure time of vessel V j should not be limited by the berthing time of vessel V k . Let M be an arbitrary large constant, then this can be modeled as for f + 1 ≤ j ≤ f + r, 1 ≤ k ≤ f + r, j k. For example, in Figure 10, it can be seen that vessels V j and V m satisfy the inequality in (31).
where E[Q(x; y; α; β, ω)] = Ω Q(x; y; α; β, ω)P(dω), and Q(x; y; α; β, ω) is the minimum value of the problem min f +r • To save fuel, the sailing speed of each vessel should be adjusted so that its arrival time at the terminal is in a given interval. This can be written as a i ≤ a i ≤ a i for i = 1, 2, . . . , f + r. • Finally, to avoid duplication and the constraints (29)-(31) among vessels, we add the constraint 1 ≤ α ik + α ki + β ik + β ki ≤ 2, for 1 ≤ i, k ≤ f + r, i < k.
We have two objectives: Minimizing the total departure delay and minimizing the fuel consumption. The total departure delay is given by where κ + := max {κ, 0} for κ ∈ R. In the model, this objective is scaled with a weight parameter, say λ, for convenience. The fuel consumption is given by [65] f i=1 c i a i +ĉ i n This function is obtained for each vessel using regression analysis. Here, c i ,ĉ i , c j (ω) andĉ j (ω) are the regression coefficients, n i (respectively, n j ) is the distance of vessel V i (respectively, vessel V j ) from the terminal, and µ i , µ j ∈ {3.5, 4, 4.5} for each i = 1, 2, . . . , f and j = f + 1, f + 2, . . . , f + r.
By introducing the auxiliary variables t i , t j , q i , q j ≥ 0, for i = 1, 2, . . . , f and j = f +1, f +2, . . . , f +r, and adding the constraints we can rewrite the summations in (34) and (35) as f +r j= f +1 c j a j (ω) +ĉ j n µ j j q j , respectively. Given the above modeling framework, the two-stage optimal random berth allocation problem is formulated as the following two-stage stochastic mixed-binary nonlinear programming model (32,33).
The nonlinear terms in (32) and (33) appear in the seventh set of constraints in each, and they can be handled by transferring them into hyperbolic constraints and then using (2) to rewrite the resulting hyperbolic constraints as a group of second-order cone constraints as follows [65]: When µ i = 3.5, the inequality 1 ≤ a µ i −1 i q i is 1 ≤ a 5/2 i q i with a i , q i > 0. Squaring both sides, we get 1 ≤ a 5 i q 2 i , which is equivalent to hyperbolic constraints min c T x + p(1)d(1) T y 1 + p(2)d(2) T y 2 + · · · + p(|Ω|)d(|Ω|) T y |Ω| s.t. Ax Using (2), the hyperbolic constraints in (37) are equivalent to the second-order cone constraints: Similarly, when µ i = 4, the inequality 1 ≤ a 3 i q i is equivalent to the following group of second-order cone constraints and when µ i = 4.5, the inequality 1 ≤ a 7/2 i q i is equivalent to the following group of second-order cone constraints Constraints (38)- (40) were obtained for each i = 1, 2, . . . , f . Similar constraints can also be obtained for each j = f + 1, f + 2, . . . , f + r. Similar to other models built in this paper, in light of the above conic transformation, a solution for Model (32,33) is now possible by applying the decomposition algorithm developed in [54].

IV. ALGORITHMS
As mentioned in Section II, there are available solution methods for SSOCP, SMILP, and DMISOCP. These optimization problems are algorithmically more mature than the more general case, the SMISOCP problem. Because the present paper captures a variety of its important applications, in this section we present a solution method for SMISOCP, leaving the development of more specific methodologies for a future research.
In the modeling process, the input requires that we solve the stochastic program for a finite number of scenarios. For this reason, in practice, the case of interest is when random parameters have a finite number of realizations. Therefore, we assume that the event space Ω is discrete, and has finitely many outcomes with known probabilities.
Problem (36) is a huge-scale DMISOCP problem, and therefore, we may not even able to solve Problem (36) using state-of-the-art optimization solvers. However, the advantage that makes the formulation in (36) different from other formulations that involve uncertainty is the nice structure that can be seen inside. We can use the solution methods available in the literature for DMISOCP, particularly cuts and relaxations, to solve the resulting extensive formulation, and this can be successfully achieved by exploiting the special structure of Problem (36). In this connection, we emphasize that all application models proposed in Section III can be formulated in the form of Problem (36).
In Subsection IV-A, we give a brief review for the solution algorithms available for solving SSOCP with more elaboration on the homogeneous self-dual method. In Subsection IV-B, we describe a solution algorithm for solving SMISOCP using cuts and relaxations by combining the existing method in Subsection IV-A for SSOCP with extensions of DMISOCP.

A. HOMOGENEOUS INTERIOR-POINT METHODS FOR SSOCP
Numerous interior-point algorithms for SSOP have been developed recently. We briefly outline these algorithms ordered chronologically. In 2014, Alzalg [30] proposed a homogeneous self-dual interior-point algorithm was proposed for SSOCP. Later in the same year, Alzalg [28] proposed a solution to the problem based on a logarithmic barrier decomposition-based interiorpoint algorithm (see also [27]). In 2015, Alzalg [29] proposed a volumetric barrier decomposition-based interior-point algorithm for solving the problem. In 2018, Alzalg et al. [31] proposed an infeasible self-dual interior-point algorithm for solving the problem. While the methodologies in [30] and [31] are based on the deterministic equivalence, which is the extensive form min c T x + d(1) T y 1 + d(2) T y 2 + · · · + d(|Ω|) T y |Ω| s.t. Ax whose dual is the problem of a stochastic program that forms an equivalent large one-stage problem containing all constraints and all scenarios, the methodologies in [28], [29] are based on Bender's decomposition, which decomposes a stochastic program into stages where, at each stage, variables at preceding stages are considered as constraints so that the subproblem at the current stage is easier to solve.
In this part, we give a brief overview of the homogeneous self-dual algorithm [30] used to solve the underlying SSOCP. The continuous relaxation of (36) is obtained by the same formulation as (36), but with x ∈ R n and y 1 , y 2 , . . . , y |Ω| ∈ R m . If we re-define d(k) to be p(k)d(k) for each k = 1, 2, . . . , |Ω|, the SSOCP is the Problem (41).
Note that we associate a different decision variable y k with each realization to be slightly more general. The dual of (41) is Problem (42).
The homogeneous model for the pair (41,42) is The following system defines the search direction system corresponding to (43).
where γ and η are two parameters, the product "•" is the Jordan multiplication defined in (1), and r p0 = p + τb − Ax, r pk = q k + τh(k) − T(k)x − W(k)y k , k = 1, . . . , |Ω|, We point out that a scaling may be needed to guarantee that we iterate on the interior of the underlying second-order cone. Common examples of such a scaling are the HRVW/KSH/M direction (introduced by Kojima et al. [66] and Helmberg et al. [67] independently and then rediscovered by Monteiro [68]), the dual HRVW/KSH/M direction (introduced by Kojima et al. [66] and rediscovered by Monteiro [68]), and the NT direction (introduced by Nesterov and Todd [69], [70]).
The work in [30] describes an efficient method to find the search directions as a solution to (44). This method exploits the special structure in (41,42) and decomposes into |Ω| smaller computations which can be performed in parallel. The results in [30] show that this method is both theoretically and computationally efficient. In particular, the number of arithmetic operations in each iteration when computing the search direction grows cubically as a function of n, m, k and l, and linearly as a function of |Ω| (see [30]).

B. GOMORY CUTS AND TIGHT RELAXATIONS FOR SMISOCP
In this part, we introduce two algorithmic developments for SMISOCP. The first development is based on the work of Cezik and Iyengar [35], which proposes Gomory cuts for solving (deterministic) mixed-integer conic programs by extending well-known techniques from mixed-integer linear programming to mixedinteger conic programs. In this part, we particularly extend the Chvátal-Gomory approach for generating cuts to the equivalent deterministic formulation (36). We introduce the following feasibility sets.
(46) The extension in this part is based on the two equivalences given in (45) (see also [35]).
Let a i and w j (k) denote the ith column of A and jth column of W(k), respectively, for k = 1, 2, . . . , |Ω|. We have the following workflow which can be applied iteratively to deterministically compute the convex hull of the feasibility set F of the formulation (see also [35,Section 2]). Workflow 1 (Chvátal-Gomory Procedure): We consider the following three steps for our setting: (i) Choose u ⪰ r 0 and v k ⪰ s 0 for k = 1, 2, . . . , |Ω|. Then (ii) Without loss of generality, assume that x ≥ 0 and y k ≥ 0, k = 1, 2, . . . , |Ω|.
The following proposition is due to [35,Lemma 1]. Proposition 1: Assume that the feasibility set F is bounded. Then every valid inequality for the convex hull of F can be obtained by repeating Workflow 1 a finite number of times.
Sherali-Adams [71], [72] and later Laserre [73] introduced hierarchies of the relaxations to present tighter relaxation of the underlying feasibility set. To extend their work to our setting, let u be a vector indexed by the empty set , subsets R ⊆ P, and sets of the form R ∪ {i}, R ⊆ P, i P. Let also v be a vector indexed by the empty set , subsets S ⊆ Q, and sets of the form S ∪ {j}, S ⊆ Q, j Q. The vectors u ∈ R (n−p+1)2p and v(k) ∈ R (m−q+1)2q are defined in (48).
For all subsets I ⊆ P, and subsets J ⊆ Q, let the values ξ I 0 , ζ I 0 (k) ∈ R, and the vectors ξ I ∈ R n , ζ J (k) ∈ R m be defined as in (48).
Finally, we apply the solution technique described in Subsection IV-A to solve Problem (49).

V. CONCLUSION
Stochastic mixed-integer second-order cone programming is an important class of optimization problems that includes stochastic second-order cone programming and stochastic mixed-integer linear programming as special cases. In this paper, we have described five different applications that lead to two-stage stochastic mixed-integer second-order cone programming, and this has been achieved by relaxing assumptions that include deterministic information to incorporate random data. Application areas include facility location, portfolio optimization, uncapacitated inventory, battery swapping stations, and berth allocation planning. It is important to indicate that all of the optimization models that have been developed in this study are intractable for an extensive form formulation. Therefore, solutions for such application models are now possible by applying either the decomposition algorithm developed recently in [54] or the scenario-based cuts proposed recently in [55]. Developing other and more specialized algorithms to solve the generic two-stage stochastic mixed-integer second-order cone programming problem is an important subject of study. In this paper, we have started this study by describing solution methods for solving the generic stochastic mixed-integer second-order cone programming problem using cuts and relaxations by combining the existing homogeneous interior-point algorithm for stochastic second-order cone programming with extensions of mixed-integer second-order cone programming. We emphasize that the work presented in this paper is far from being a survey or a comprehensive study, but we believe that it would rather shape a basis for further research studies in secondorder cone programs over integers under uncertainty. In particular, it is our firm belief that this research is a rich source for stochastic mixed-integer second-order cone programming models, and it can be greatly utilized to implement numerous future algorithms for solving this class of optimization problems.