Aircraft Loading Optimization: MemComputing the 5th Airbus Problem
AAircraft Loading Optimization: MemComputing the th Airbus Problem
Fabio L. Traversa ∗ MemComputing, Inc., San Diego, CA, 92130 CA (Dated: June 12, 2019)On the January 22nd 2019, Airbus launched a quantum computing challenge to solve a set ofproblems relevant for the aircraft life cycle (Airbus challenge web-page). The challenge consistsof a set of 5 problems that ranges from design to deployment of aircraft. This work addressesthe 5 th problem. The formulation exploits an Integer programming framework with a linear ob-jective function and the solution relies on the MemComputing paradigm. It is discussed how touse MemCPU TM software to solve efficiently the proposed problem and assess scaling properties,which turns out to be polynomial for meaningful solutions of the problem at hand. Also discussedare possible formulations of the problem utilizing non-linear objective functions, allowing for differ-ent optimization schemes implementable in modified MemCPU software, potentially useful for fieldoperation purposes. I. INTRODUCTION
Automated computation in the 21st century hasreached a paramount role in industry, finance, consumertechnology and much more [1, 2]. However, the morethat computation becomes relevant, the greater the chal-lenges presented to both industry and academia. So-lutions to today’s computational limits are provided bymore and more sophisticated and high-performing CPUs,GPUs etc. [3], as well as by the insurgence of distributedcomputing in cloud infrastructures [4].Nevertheless, the unceasing growth of computing de-mand is pushing towards completely new architecturesas well as new paradigms. This can mean not only “han-dling more computation,” but also making the computa-tion more efficient. For example, neuromorphic comput-ing based on diverse realizations of artificial neural net-works [5–7] promises a paradigm shift in machine learningand artificial intelligence.In this scenario, new computing architectures basedon quantum physics and not strictly on algorithmic ap-proaches promise to solve among the most challengingproblems ranging from drug discovery to A I[8–11]. How-ever, it is still not clear if and when reliable quantumcomputers will show what currently goes under the nameof quantum supremacy over classical computing architec-tures [12–15].However, companies and academia are already lookingat applications for future quantum computers [13, 16, 17].I focus here on the challenge recently launched by Air-bus [17] consisting of 5 problems related to the life cycleof an aircraft, i.e., ranging from the design to the de-ployment of an aircraft. The general requirements forthe challenge are that for each problem a suitable algo-rithm implementable on a quantum computer should bedeveloped, tested on a simulated quantum computer, andfinally assessed in terms of the performance at scale.I consider in this work the 5th Airbus problem [17]:the Aircraft Loading Optimization (ALO) problem. The ∗ email: [email protected] ALO requires optimizing the placement of containers ofdifferent sizes and weights in an aircraft subject to limi-tations on maximum weight allowed, maximum tolerableshear and center of gravity. However, instead of develop-ing on algorithm for future quantum computers, I presenta formalization of the problem using Integer Program-ming (IP) framework [18]. IP is a problem formulationwidely used in industry and academia consisting of an ob-jective function to be minimized over a set of variables.Variables may be constrained to take binary, integer orcontinuous values. Moreover, the objective function maybe subject to further constraints in the form of linearinequalities among variables. Finally, depending on thenature of the objective function, we can have differentIP problems such as linear, quadratic or other non-lineartype.The formulation I consider here covers all requirementsof the ALO statement and leads to an IP problem in-volving binary variables only and having two objectivefunctions to be minimized, one linear and the other non-linear. I therefore propose two different solution strate-gies to solve this multi-objective IP problem.The binary IP (BP) problem is NP-complete, some-times also known as , and is one ofKarp’s 21 NP-complete problems [19]. Therefore, it is notsurprising that Airbus poses a challenge about a prob-lem whose natural formulation leads to a special case ofon NP-complete problem. Several general-purpose opensource [20–24] and commercial solvers [25–28] have beendeveloped to solve IP problems. However, these can failwhen the problems are particularly hard and thereforespecialized solvers optimized to exploit specific structuresfor ILP have also been developed [29–32].IP problems can be approached employing differenttechniques. A classification of those includes heuris-tics [31, 33–35] and exhaustive algorithms [18, 36–38] alsocalled “complete algorithms.” The complete algorithmfor ILP that is most commonly employed is a combina-tion of cutting planes and branch-and-bound, also knownas the branch-and-cut algorithm [18]. All these solvershave demonstrated varied degrees of success on a varietyof IP problems [39–43], however scaling properties for BP a r X i v : . [ c s . ET ] J un xL/2-L/2 01 2 3 …. …. NN-1N-2
12 1 cgmax x cg min x cg F F xL/2-L/2 0 S max (x) S(x)
FIG. 1: Sketch of the Airbus Aircraft Loading Optimization problem. The aircraft is to be loaded by selecting aportion of the available payload. The containers can have different shapes. For this ALO, 3 different shapes areconsidered numbered from 1 to 3 as in the figure. The maximum weight allowed is limited. The center of gravity x cg of the system aircraft + containers must remain within the limit [ x min cg , x max cg ] and the shear curve S ( x ) must also belimited to remain under a give shear limit curve S max ( x ).problems such as the one proposed by Airbus can easilyshow exponential blowing-up when applying one of thosemethods [19].In this work I discuss the general purpose approachto the solution of the Aibus ALO problem based on the memcomputing paradigm introduced in [44]. Memcom-puting (computing with and in memory [45]) approachis neither based on stochastic search nor on trial and er-ror strategy. In addition, memcomputing does not use aset of instructions recursively employed to find solutionsto the problem at hand. Therefore, we can classify thememcomputing approach as non-algorithmic [45].The memcomputing approach relies on embedding agiven problem into an electronic circuit, which repre-sents a possible realization of a memcomputing machine(MM) [44–46]. These electronic circuits, when designedin order to satisfy crucial properties [46], naturally re-lax to an equilibrium that represents the solution of theproblems at hand [45–48].In order to approach the Airbus ALO problem I usethe MM realization described in [35] employing self-organizing algebraic gates (SOAGs). The SOAG-basedMM can be efficiently simulated in software and in thiswork I use the software developed at MemComputingInc. [49] named MemCPU coprocessor.In section II I briefly describe the Airbus ALO prob- lem. In section III the BP formulation of the ALO prob-lems is discussed. In section IV a quick overview of mem-computing and MemCPU software is provided, while insection V are discussed both numerical results and scal-ing properties of MemCPU software applied to the ALOproblem. II. AIRBUS AIRCRAFT LOADINGOPTIMIZATION PROBLEM
The ALO problem proposed by Airbus is sketched inFig. 1. Following the problem statement we need to opti-mize the placement of a subset of containers picked froman available payload formed by n different containers.Each container is described by the triplet ( k, m k , s k ) with k the identification number, m k the mass and s k the sizeof the k − th container.The aircraft has N available positions for standardcargo as depicted in Fig. 1. The size can be of 3 dif-ferent types: size 1 occupies one position; size 2 occupieshalf a position and may share the position with anothersize 2 container; size 3 occupies 2 positions.The problem requires maximizing the mass of the car-ried freight for the flight such that:(A) the containers must be placed consistently with sizeand positions available(B) the mass of the carried freight does not exceed themaximum payload capacity of the aircraft W p (C) the center of gravity position of the carried freight+ aircraft must be within the limits [ x min cg , x max cg ](D) the shear curve S ( x ) defined by the container distri-bution must be bounded by a given maximal shearcurve S max ( x )(E) while maximizing the mass, the position of the centerof gravity of the carried freight + aircraft should beoptimized to be as close as possible to a target centerof gravity x tcg .These constraints are sketched in Fig. 1. While theconstraints (A) and (B) do not need additional details orexplanations, the others do. In (C) the center of grav-ity is defined by the the equilibrium of all forces (con-tainer + aircraft) in the gravity axis direction (see force F in Fig. 1). Each container is supposed to have uniformmass, therefore the weight force can be considered as theweight of a point located at the center of the container.Finally, the mass and location of the center of gravity ofthe empty aircraft are provided.The constraint (D) requires a consistent definition ofthe shear curve. The shear curve is defined by S ( x ) = (cid:90) x − L m ( x (cid:48) ) dx (cid:48) for x < a ) S ( x ) = (cid:90) L x m ( x (cid:48) ) dx (cid:48) for x > b )where L is the length of the loading region of the aircraftand x is the position relative to the center of the loadingzone (see Fig. 1). In the statement of the problem, themaximum shear curve S max ( x ) can be either linear andsymmetric with respect to x = 0, or asymmetric or somenon-linear function of x .Finally, (E) represents an extra optimization require-ment and not really a constraint. In fact, from the Air-bus statement, what should be primarily optimized isthe mass of the carried freight, then optimize the distri-bution such that the resulting center of gravity is as closeas possible to a target. III. INTEGER PROGRAMMINGFORMULATION
A natural way to mathematically formulate the prob-lem posed in section II is using Integer Programmingframework [18].Let us start assigning a binary variable y k,j to eachcontainer k = 1 , ..., n located in the position j = 1 , ..., N of the aircraft. Since each container can be located in at most one position in the aircraft we constrain thesevariables requiring N (cid:88) j =1 y k,j ≤ k = 1 , ..., n. (2 a )On the other hand, we have the constraint on the sizeand the number of available spots on the aircraft thatcan be expressed as a non-linear inequality. In order toformalize this constraint, let us further characterize thevariable y k,j . We consider here, for simplicity, only the3 sizes defined in the Airbus statement. However, thisformulation can be easily extended to more complicatedsizes and bin distributions in the aircraft. Let us define K , K and K as three sets of indexes such that K ∪ K ∪ K = { , ..., n } . Moreover, if k ∈ K h then s k = h .Therefore, K , K and K regroup the variables in setsof variables corresponding to the same container size.Using this index characterization, the constraint (A)on the sizes and bins can be expressed as (cid:88) k ∈ K y k, + 12 (cid:88) k ∈ K y k, + (cid:88) k ∈ K y k, ≤ b (cid:48) ) (cid:88) k ∈ K y k,j + 12 (cid:88) k ∈ K y k,j + (cid:88) k ∈ K ( y k,j − + y k,j ) ≤ j = 2 , ..., N − b (cid:48)(cid:48) ) (cid:88) k ∈ K y k,N + 12 (cid:88) k ∈ K y k,N + (cid:88) k ∈ K y k,N − ≤ b (cid:48)(cid:48)(cid:48) )It is easy to prove that these inequalities guaran-tee that for each bin there is no overlapping of con-tainers except for containers of size 2 for which two ofthem can occupy the same bin (the coefficient 1 / (cid:80) k ∈ K ( y k,j − + y k,j ) in which two consecutive containersof size 3 appear in the same constraint; therefore a con-tainer of size 3 occupies 2 bins if selected. It is also worthnoticing that containers of size 3 have j that ranges onlyfrom 1 to N − N − (cid:88) k,j m k y k,j ≤ W p . (2 c )The constraint (C) on the center of gravity requiresthe definition of signed distance d s k ,j from x = 0 for acontainer of size s k located in the bin j . Consideringequally distributed bins around x = 0 (this assumptionis the same reported in the Airbus statement but can beeasily relaxed), we have that d s k ,j = 2 j − N − N L for s k = 1 , a ) d s k ,j = 2 j − N N L for s k = 3 (3 b )Using this distance definition, the center of gravity ofthe loaded aircraft can be evaluated as x cg = (cid:80) k,j m k d s k ,j y k,j + W e x ecg (cid:80) k,j m k y k,j + W e (4)where W e and x ecg are respectively the mass and the cen-ter of gravity of the empty aircraft.Eq. (4) allows the formalization of the constraint (C)through the inequalities (cid:88) k,j m k ( d s k ,j − x max cg ) y k,j ≤ W e ( x max cg − x ecg ) (2 d (cid:48) ) (cid:88) k,j m k ( x min cg − d s k ,j ) y k,j ≤ W e ( x ecg − x min cg ) (2 d (cid:48)(cid:48) )that express x cg ≤ x max cg and x cg ≥ x min cg respectively.Regarding the constraint (D), the shear curve can beeasily calculated at each bin location and set smaller than S max ( x ): (cid:88) k ∈ K ∪ K j (cid:88) j (cid:48) =1 m k y k,j (cid:48) ++ (cid:88) k ∈ K j − (cid:88) j (cid:48) =1 m k y k,j (cid:48) + 12 m k y k,j ≤≤ S max ( x j ) for j ≤ N/ e (cid:48) ) (cid:88) k ∈ K ∪ K N (cid:88) j (cid:48) = j m k y k,j (cid:48) ++ (cid:88) k ∈ K N − (cid:88) j (cid:48) = j − m k y k,j (cid:48) + 12 m k y k,j ≤≤ S max ( x j ) for j ≥ N/ e (cid:48)(cid:48) )which expresses the shear curve defined in Eq. (1) at thecenters of the bins x j taking into account the differentsizes of the containers.Finally we can define the objective function of the ALOproblem as f ( y ) = − (cid:88) k,j m k y k,j . (5)whose minimization over y , subject to the collection ofconstraints (2), provides the distribution and the maxi-mum mass of the carried freight satisfying the constraints(A)-(D) of section II. A. Center of Gravity Optimization
We briefly discuss in this section the further opti-mization required in (E) of section II. This can be han-dled with the following scheme. Once the problem (5)is solved, we obtain the maximum mass of the carriedfreight as W max = (cid:88) k,j m k ˜ y k,j (6)where ˜ y is the solution of the minimization of (5) subjectto (2). We therefore define a new IP problem where weinclude the constraints (2a,2b,2c,2e) but exclude for now(2d). Instead we include the extra constraint − (cid:88) k,j m k y k,j ≤ − τ W max . (7)where 0 ≤ τ ≤ is a tolerance parameter that can beused in case we are interested in a slightly lighter carriedfreight but better center of gravity.We then define the new objective function f ( y ) = (cid:26) x cg − x tcg if x cg ≥ x tcg − x cg + x tcg if x cg ≤ x tcg (8)with x cg defined by Eq. (4).The target now is to find the assignment of y thatminimizes Eq. (8) under the constraints (2a,2b,2c,2e) and(7). This can be achieved either using a sequence of linearIP problems in which we reintroduce the constraint (2d)and we iteratively shrink the [ x min cg , x max cg ] around x tcg , orby implementing it directly as a non-linear IP problemwith objective function Eq. (8). IV. MEMCOMPUTING APPROACH
The memcomputing approach to IP problems isbased on the concept of Self-Organizing Algebraic Gates(SOAGs) introduced in [35]. SOAG is a novel circuitdesign developed at MemComputing, Inc. [49] by the au-thor of this work and is part of a class of self-organizingcircuits as Self-Organizing Logic Gates (SOLGs) intro-duced in [46, 50].Both SOLGs and SOAGs are building blocks for prac-tical realizations of Universal Memcomputing Machines(UMM) [44, 45, 51], with digital input-output (DigitalMM, DMM) [46].The SOAG is designed to self-organize toward an al-gebraic relation representing a linear inequality amongvariables. In this work, the SOAGs are designed to sat-isfy linear relations between binary variables (see Fig. 2)since the ALO problem formulation involves binary vari-ables only.By connecting together SOAGs we form a Self-Organizing Algebraic Circuit (SOAC), see Fig. 3. The
Self-Organizing Algebraic Gate 𝑗 𝑎 𝑗 𝑥 𝑗 − 𝑏 𝑗 𝑎 𝑎 𝑎 𝑎 𝑜𝑢𝑡 DCM DCMDCMDCM DCM
Dynamic correction module
DCM
Read Voltages
Inject Current
FIG. 2: Reprint from [35]. Sketch of a Self-OrganizingAlgebraic Gate. All terminals allow a superposition ofincoming and outgoing signals from the surrounding cir-cuit. The central unit processes the signals in order tosatisfy a linear algebraic relation consistent with the re-quirement of the “out” terminal. The self-organization isenforced by the Dynamic Correction Modules that readvoltages from all terminals and inject a current to theappropriate terminal as long as the algebraic relation isnot satisfied. 𝑥 𝑥 𝑥 𝑥 𝑠𝑎𝑡𝑖𝑠𝑓𝑦 𝑠𝑎𝑡𝑖𝑠𝑓𝑦 𝑠𝑎𝑡𝑖𝑠𝑓𝑦𝑠𝑎𝑡𝑖𝑠𝑓𝑦 𝑠𝑎𝑡𝑖𝑠𝑓𝑦 𝑠𝑎𝑡𝑖𝑠𝑓𝑦 𝑥 𝑥 𝐴 𝑒𝑞 𝑥 = 𝑏 𝑒𝑞 𝐴 𝑖𝑛𝑒𝑞 𝑥 ≤ 𝑏 𝑖𝑛𝑒𝑞 FIG. 3: Reprint from [35]. Sketch of a Self-OrganizingAlgebraic Circuit (SOAC). SOAGs are connected to-gether in an architecture that directly maps the IP intothe SOAC.SOAC collectively self-organizes in order to satisfy therelations embedded in the gates. In this way, it is trivialto embed an IP problem directly into the SOAC. Eachinequality in (2) or (7) can be mapped directly into aSOAG. The objective functions (5) can be easily refor-mulated as an extra linear inequality (cid:88) k,j f k,j y k,j ≤ ˜ b (9)where ˜ b represents a threshold for the maximum carriedfreight that is dynamically changed each time a feasi-ble solution is found in order to find solutions closerand closer to the global minimum of the problem. On the other hand, the objective function (8) can be imple-mented as a pair of linear inequalities of the form (cid:88) j f + k,j (˜ b ) y k,j ≤ g + (˜ b ) (10) (cid:88) j f − k,j (˜ b ) y k,j ≤ g − (˜ b ) (11)where ˜ b is an extra parameter that again can be dynami-cally changed to find solutions of increasing quality, eachtime closer to the global optimum. f ± and g ± are linearfunctions of ˜ b defined by substituting (4) in (8). Thisalso leads to defining ˜ b as the threshold of the distanceof the actual center of gravity from the target x tcg .The problem formulated and embedded in a SOAC asdescribed can be efficiently handled by either actuallybuilding the electronic circuit or just by simulating itsince it involves only standard (non-quantum) electroniccomponents. Simulating means solving differential equa-tions of the form dQ ( i, v, x ) dt = F ( i, v, x ) , (12)with appropriate initial conditions and where Q and F are non-linear functions of the voltages ( v ), currents ( i )and extra internal state variables ( x ) characterizing thecircuit. In this picture a configuration of the voltages( v ) at a given time represents an actual assignment tothe variables of the IP problem by reading the voltagesthrough thresholds: if a voltage at a node of the circuitis above the threshold, then it corresponds to a logical 1,and otherwise corresponds to a logical 0. On the otherhand, the transition function of these machines (namelythe function that maps input to output) is physical (ana-log) and takes full advantage of the collective state ofthe system to process information [46, 51, 52]. However,despite its analog nature, the mapping of voltages intobinary variables through thresholds allows efficient sizescaling for these machines, avoiding the bottleneck re-lated to the precision of writing and reading inputs andoutputs. Finally, it is worth noticing here that the mem-computing approach used to solve binary problems aspresented in [35] for the case of IP problems, does notprovide proof of optimality for a given solution, nor doesit detect the infeasibility of a problem.In this work we consider the simulation of the SOACimplemented in the MemCPU software developed atMemComputing, Inc.[49] and available also as a Softwareas a Service offering.The working principle of SOAGs and SOACs, i.e., self-organization of voltages and currents of an electronic cir-cuit in order to satisfy algebraic relations, is enforcedby both active and passive electronic elements with andwithout memory [35, 45, 46]. The features of DMMs havebeen investigated, and interesting properties emerge froma correct design such as long-range order correlations andtopological robustness [53, 54]. If mathematical require-FIG. 4: Histograms of container mass from Airbus dataset reported in Table I.ments described in [46] are fulfilled, then persistent os-cillations or chaos are avoided [47, 48]. Self-organizingdigital circuits (i.e. both SOLCs and SOACs), being arealization of DMMs can in principle solve efficiently (i.e.with polynomial resources) complex combinatorial opti-mization problems, given that mathematical features onthe constitutive equation (12) are satisfied [44–46]. Self-organizing digital circuits (i.e. both SOLCs and SOACs)have also been proved to efficiently handle a variety ofcombinatorial optimization problems ranging from max-imum satisfiability (MAXSAT) [55, 56] to quadratic un-constrained binary optimization (QUBO) [57] and fromIP [35] to the training of neural networks [58]. V. SCALING RESULTS
In order to assess the scaling properties of solving theAirbus ALO problem employing MemCPU coprocessorwe need to generate a set of meaningful benchmarks atdifferent values of N and n . Since we have no informa-tion about the actual values of N and n and not even onthe ratio between number of containers at different sizesor typical mass distribution of the containers available,we use the sample data set provided by Airbus (see Ta-ble I for the data set, or it can also be found from [17])to extrapolate. Notice that the data set provided fromAirbus is just “for illustration and for testing the algo-rithm” [17].In order to generate benchmarks for different sizes weneed to generate a distribution of containers. One of thepossible ways is to try to produce a distribution similarto that from the sample data set. In Fig. 4 the distri-butions of container masses from the data set in Table Iis reported. Even though there are not enough data todetermine with certainty the distribution shape, it is rea-sonable to assume that both distributions can be recre-ated from a bimodal Gaussian distribution, cut off atcertain boundaries. N = 20 n = 30 W p = 40000 W e = 120000 x ecg = − . Lx min cg = − . Lx max cg = 0 . Lx tcg = 0 . LS max ( x = 0) = 22000, linear, symmetric k s k m k (kg)1 1 21342 1 34553 1 18664 1 16995 1 35006 1 33327 1 25788 1 23159 1 188810 1 178611 1 327712 1 298713 1 253414 1 211115 1 260716 1 156617 1 176518 1 194619 1 173220 1 164121 2 180022 2 98623 2 87324 2 176425 2 123926 2 148727 2 76928 2 83629 2 65930 2 765TABLE I: Airbus ALO data set [17].In Fig. 11 I have included Matlab code that generatesALO problems for any n and N (it is also available aconverter from Matlab to .mps format at [59]). The dis-tribution of containers is generated in the subfunction container distribution generator . A bimodalGaussian distribution of containers is generated for eachsize type. The ranges are chosen in such a way that thedistributions in Fig. 4 can be qualitatively reproduced.It is also worth noticing that the masses are scaled de-pending on N . In fact, maintaining the same weight ofthe airplane and the same limit on the maximum weight, ob j ec ti v e Time Const. 1 Time Const. 2 Time Const. 3 Time Const. 4 ob j ec ti v e Time Const. 5 Time Const. 6 Time Const. 7 Static Const. 1 ob j ec ti v e Static Const. 2 Static Const. 3 Static Const. 4 Static Const. 5 ob j ec ti v e Static Const. 6 Static Const. 7 Static Const. 8 Static Const. 10
FIG. 5: Monte Carlo distribution of the objective (5) versus MemCPU coprocessor parameters [49]. The ALO problemis defined with the parameters in Table I and x min cg = − . A. MemCPU software results
As discussed in section IV, a MemCPU software is anemulator of an electronic circuit with a given IP prob-lem embedded within. However, in order to be efficientin solving the problem at hand, the MemCPU coproces- sor needs a proper set of parameters characterizing elec-tronic elements such as resistors, capacitors, etc. Theseparameters do not depend on the size of the problem(they are scale free) but depend on the structure of theproblem. The MemComputing SaaS currently providesa Monte Carlo routine that evaluates the distribution ofthe objective function of the IP problem as a functionof the parameters [60]; however future releases will pro-vide a predictor routine for parameters that will avoidthe Monte Carlo-based tuning [60]. The objective func-tion distribution is reported in Fig. 5. This distributionhas been computed for the problem generated using thecode in Fig. 11 using the container distribution in Ta-ble I as input and x min cg = − . x min cg isonly for tuning purposes because this restricts the possi-ble feasible solutions of the problem making it “harder.”In fact, it is easy to estimate the limit of the center ofgravity x lim cg for any selection of containers ( τ = 0) fromthe available payload. In Fig. 6, the forces in the cen-ter of gravity evaluation (Eq. 4) are summarized for thelimit of containers of equal weight corresponding to themaximum weight allowed by the maximum shear curve.The value is x lim cg = − / ≈ − . x min cg = − .
006 that largely reduces the number of feasiblesolutions. This choice for x min cg leads to minimal feasiblesolutions that have objective far from maximum allowedcarried freight W p , however, used for tuning is a goodchoice because largely reduces the number of feasible so-lutions. This choice helps to produce a sharper distribu- x L/2 -L/2 max (x)=S(x)22.000 F r L/418.000 -9L/40 F l -9L/80 S(x)S max (x) 44.000/L80.000/L W e x cg -L/20 lim S max (x)=S(x) F r L/4
S(x) S max (x) W e x cg -L/20 lim L284=-̶̶ τ = 0τ = 1
FIG. 6: Configuration of forces to evaluate how close thecenter of gravity x lim cg can approach the target x tcg = 0 . τ = 0) and with maximum selectedpayload ( τ = 1). For the maximum selected pay-load, the maximum density of 80000 /L (container massper unit of aircraft length) is assumed, consistent withthe data set in Table I and the subfunction con-tainer distribution generator in Fig. 11.tion of objective function values versus parameters whichalso is useful for solving the problem in section III A.From the objective distribution of Fig.5, a set of theparameters for the MemCPU coprocessor can be easilyextracted and used for running any other instance of theALO problem. However, extracting a good set of param-eters from the distribution can be done in multiple ways.In this work, the following procedure has been followed:from each Markov chain the set of parameters returningthe best objective is selected. If there is more than oneset of parameters with the same objective for the sameMarkov chain, the choice is made based upon the onethat has the best average of the objectives related to theclosest 10 sets of parameters within the same Markovchain. This super-selection rule comes from the fact thateach MemCPU coprocessor run starts with random ini-tial conditions for voltages and currents, and thereforethere could be a fluctuation in the output objective dueto the finite (and short) simulation time [60]. Once wehave the best parameter choice per Markov chain, we canrun the same problem multiple times with random ini-tial conditions for the same set of parameters and selectthe set of parameters that statistically arrives fastest to
20 30 40 50 60 70 80 90 100 N -1 ti m e , s FIG. 7: MemCPU coprocessor run time versus N to finda configuration of containers with total mass larger thatthan 99.9% of the maximum mass allowed on the air-craft or of the total payload if smaller than the maxi-mum mass allowed. Dots connected by solid lines are theaverage time for 100 different ALO instances generatedusing the code in Fig. 11. Dashed curves are the scalingrelation (13). Runs have been carried out on an Intel(R)Xeon(R) Gold 6138 CPU @ 2.00GHz. n l -1 ti m e , s FIG. 8: MemCPU coprocessor run time versus n l to finda configuration of containers with total mass larger thatthan 99.9% of the maximum mass allowed on the air-craft or of the total payload if smaller than the maxi-mum mass allowed. Dots connected by solid lines are theaverage time for 100 different ALO instances generatedusing the code in Fig. 11. Dashed curves are the scalingrelation (13). Runs have been carried out on an Intel(R)Xeon(R) Gold 6138 CPU @ 2.00GHz.the best objective in a given time out. In this way, Iextracted the set of parameters reported in Table II.Using the parameters in Table II, the MemCPU soft-ware has been tested on a set of ALO instances gener-ated for different n and N using the code in Fig. 11.The created benchmark consists of 100 ALO instancesfor each pair ( n, N ). Each instance has a differentcontainer mass distribution randomly generated by thesubfunction container distribution generator of Fig. 11. For each instance, the container sizes were n = dim( K ) = n/ n = dim( K ) = n/ n = dim( K ) = n/ n j and n + n + n = n .Since the value of the objective related to the globalminimum of the ALO problem is not known and theMemCPU technology does not provide proof of optimal-ity, I have set a very tight threshold to assess the qualityof the solution: accept the solution if the objective func-tion is smaller than − τ W p if W p ≤ (cid:80) m k or smaller than − τ (cid:80) m k otherwise, with τ = 0 . τ = 0 .
999 is reported as a function of N (Fig. 7) or as afunction of the number n l of nonzero entries in the ma-trix defined by constraints (2) (Fig. 8) for different ratios r = n/N .The results as a function of n l are useful for properlyassessing the scaling of the MemCPU coprocessor. In-deed, the MemCPU software associates an SOAG with anumber of terminals equal to the terms in the constraintto each constraint. Each terminal includes a simulateddynamic correcting module (DCM) [35, 46]. Therefore,the number of DCMs that the MemCPU software sim-ulates is equal to n l , and it is natural to measure thecomplexity in terms of n l that is at the same time also ameasure of the size of an IP problem. The interpolationcurve on log-log scale in the plane n l -time provides thefollowing relation: t = 10 − . r − . n . r +1 . l (13)where the dependence on r is valid for 0 . ≤ r ≤
3. Theequation (13) shows a sub-quadratic scaling of the timeto solution versus n l . From (2) it is easy to realize that n l ∝ nN and, substituting this relation in equation (13),the scaling as a function of N and n can be recovered.It is worth discussing the dependence of Eq. (13) on r . The first thing that can be noticed is that the de-pendence of n l on r is very weak. The coefficient 0.11in (13) makes the exponent range from 1.26 to 1.58 for0 . ≤ r ≤
3. On top of the weak dependence, there iscompensation from the stronger and negative dependenceof the prefactor 10 − . r − . on r . However, increasing r ,an evident dependence of the exponent on r disappears,and it saturates for r >
3. This is not surprising behaviorbecause for larger r the problem becomes “easier” sincethere are many more choices of container selections thatrespect the constraints, and the burden in the calcula-tion depends only on the size of the ALO. We do notshow here numerical experiments to support the latterclaim because considering n much larger the N seemsmeaningless in practical cases.
20 30 40 50 60 70 80 90 100 N -1 ti m e , s FIG. 9: MemCPU coprocessor run time versus N to finda configuration of containers that a) has a total masslarger that than 99.8% of the maximum mass allowed onthe aircraft or of the total payload if smaller than themaximum mass allowed and b) optimizes the center ofgravity of the aircraft. Dots connected by solid lines arethe average time for 100 different ALO instances gener-ated using the code in Fig. 11. Dashed curves are thescaling relation (13). Runs have been carried out on anIntel(R) Xeon(R) Gold 6138 CPU @ 2.00GHz.Let us consider now the ALO problem described in(section III A). In order to handle this problem by usinga non-modified version of MemCPU software, I have cre-ated a sequence of linear IP problems that converges tothe solution of the non linear problem of (section III A).This sequence of problems can be easily created makinguse of the code of Fig. 11 as follows:a. Set x min cg and x max cg for some initial values, for example,the ones given in Table I.b. Generate the problem for a given x min cg and x max cg .c. Add to the matrix of the constraints the objectivefunction of the problem as an extra constraint. For theright hand side of this constraint set − τ W max where W max = W p if W p < (cid:80) m k or W max = (cid:80) m k other-wise.d. Substitute the objective function with a null objectivefunction and solve the problem.e. If a solution to the problem is found, evaluate x cg bymeans of J , VL , V , xe cg and We as reported in thecode of figure 11f. if x cg < x tcg set x min cg = x cg + (cid:15) otherwise set x max cg = x cg − (cid:15) for some (cid:15) >
0. Go back to b.In Fig. 9 the time to solution is reported for selectingthe payload that maximizes the mass ( tau = 0 .
998 hasbeen used in this case) and at the same time providesthe center of gravity as close as possible to the target,0
20 30 40 50 60 70 80 90 100 N -0.06-0.04-0.0200.020.040.060.080.10.12 0.5 0.6 0.7 0.8 0.9 n/N m a ss , kg FIG. 10: Center of gravity of the system aircraft + se-lected payload with total carried mass larger than 99 . W max = min( W p , (cid:80) m k ). In the inset, the meanof the total mass generated by the subfunction con-tainer distribution generator of Fig. 11 and W p versus r = n/N is plotted.for different values of N and n . The dashed curves thatinterpolate this version of the ALO result as a shift ofthe equation (13) where the coefficient 10 − . r − . is re-placed by 10 − . r − . . Finding the same scaling proper-ties as in the other problem is not surprising, since solv-ing this ALO problem follows the steps a.-f., and so is acascade of problems similar to the ones solved in Fig. 7.Finally, it is worth noticing where the optimized cen-ters of gravity are located. In Fig. 10, we can see that for n = N/ W max ≈ W p / n ≥ N , W max > W p (see inset of Fig. 10) and the optimized centers of gravityapproach the second limit of Fig. 6. B. Further improvements
We briefly discuss here a few improvements that canbe made in order to have an even more efficient time-to-solution for the Airbus ALO problem computed througha MemCPU coprocessor. However, since they do notqualitatively affect the scaling assessed in this work, theimplementation of these is beyond the current scope.The first improvement concerns the initial values ofthe objective function (5). The MemCPU coprocessorsearches solutions that have an objective equal to or lowerthan a target. Once found it decreases the objective andsearches for the next solution (see section IV and [35]).In this work I have used infinite as the initial objective. With this choice, the MemCPU software provides manysolutions at lower and lower objectives, progressing fromthe larger to smaller values. Even though this conver-gence is usually fast, this can be improved even moreby giving as a starting value for the objective function τ min( W p , (cid:80) m k ).Another improvement that scales the time to solutiondown by a few orders of magnitude is running the Mem-CPU coprocessor on GPUs rather than CPUs [61]. Infact, the computation performed by MemCPU coproces-sor is nothing other than a circuit simulation that canefficiently be distributed on GPUs [61].Finally, the search for maximal carried freight that op-timizes the center of gravity at the same time can beaccelerated by directly implementing the nonlinear ob-jective function (8), avoiding the sequence of ALO prob-lems discussed in section V A and solving only one ALOproblem. VI. CONCLUSIONS
In summary, I have shown how to employ MemCom-puting through the the MemCPU coprocessor to tacklethe 5th Airbus problem efficiently without a quantumcomputer.In order to use the the MemCPU coprocessor for the5th Airbus problem, which is an aircraft loading opti-mization problem, an IP problem has been formulatedthat includes all constraints required by Airbus with noexception and no approximation.The scaling properties assessed in this work show sub-quadratic scaling by the MemCPU coprocessor as a func-tion of the size of the problem measured by the numberof nonzero elements of the constraint matrix of the asso-ciated IP problem.The scaling properties of the MemCPU coprocessor al-low efficient solutions of large ALO problems and it rep-resents a solution ready to be deployed in the field.Finally, I have also discussed extra improvements thatcan be made in order to further (and substantially) ac-celerate the time to solution for the ALO problem by theMemCPU coprocessor.
VII. ACKNOWLEDGMENTS
I thank the MemComputing team for useful discussionsand in particular Dr. Tristan Sharp and John Beanefor the revision of the manuscript. This work has beensupported by MemComputing, Inc.
VIII. REFERENCES [1] Eiben, A. & Smith, J. Introduction to Evolutionary Com-puting (Springer Berlin Heidelberg, 2015).[2] Kephart, J. & Chess, D. The vision of autonomic com-puting.
Computer , 41–50 (2003).[3] Vestias, M. & Neto, H. Trends of CPU, GPU and FPGAfor high-performance computing. In (IEEE, 2014).[4] Hashem, I. A. T. et al. The rise of “big data” on cloudcomputing: Review and open research issues.
Informa-tion Systems , 98–115 (2015).[5] Burr, G. W. et al. Neuromorphic computing using non-volatile memory.
Advances in Physics: X , 89–124(2016).[6] Torrejon, J. et al. Neuromorphic computing withnanoscale spintronic oscillators.
Nature , 428–431(2017).[7] Furber, S. B., Galluppi, F., Temple, S. & Plana, L. A.The SpiNNaker project.
Proceedings of the IEEE ,652–665 (2014).[8] Linke, N. M. et al.
Experimental comparison of two quan-tum computing architectures.
Proceedings of the NationalAcademy of Sciences , 3305–3310 (2017).[9] Bennett, C. H. & DiVincenzo, D. P. Quantum informa-tion and computation.
Nature , 247–255 (2000).[10] Ladd, T. D. et al.
Quantum computers.
Nature ,45–53 (2010).[11] Raha, K. et al.
The role of quantum mechanics instructure-based drug design.
Drug Discovery Today ,725–731 (2007).[12] Harrow, A. W. & Montanaro, A. Quantum computa-tional supremacy. Nature , 203–209 (2017).[13] MIT Technology Review.[14] Dyakonov, M. When will useful quantum computers beconstructed? not in the foreseeable future, this physicistargues. here’s why: The case against: Quantum comput-ing.
IEEE Spectrum , 24–29 (2019).[15] Tang, E. A quantum-inspired classical algorithm forrecommendation systems. arXiv:1807.04271 (2018).http://arxiv.org/abs/1807.04271v1.[16] Accenture report.[17] Airbus challenge press-release; Airbus challenge web-page.[18] Schrijver. Theory of Linear Integer Programming (JohnWiley & Sons, 1998).[19] Garey, M. R. & Johnson, D. S.
Computers and In-tractability; A Guide to the Theory of NP-Completeness (W. H. Freeman & Co., New York, NY, USA, 1990).[20] Gnu linear programming kit, version 4.32. URL .[21] Forrest, J. et al.
Coin-or/cbc: Version 2.9.9 (2018).[22] Gleixner, A. et al.
The SCIP OptimizationSuite 6.0. ZIB-Report 18-26, Zuse Institute Berlin(2018). URL http://nbn-resolving.de/urn:nbn:de:0297-zib-69361 .[23] Ralphs, T. et al.
Coin-or/symphony: Version 5.6.16(2017).[24] Berkelaar, M., Eikland, K. & Notebaert, P. lp solve 5.5,open source (mixed-integer) linear programming system.Software (2004). URL http://lpsolve.sourceforge.net/5.5/ . [25] URL .[26] URL .[27] URL .[28] URL http://mathworks.com .[29] Applegate, D., Bixby, R., Chvatal, V. & Cook, W. Con-corde tsp solver (2006).[30] Floudas, C. A. & Lin, X. Mixed integer linear program-ming in process scheduling: Modeling, algorithms, andapplications.
Annals of Operations Research , 131–162 (2005).[31] Boston, K. & Bettinger, P. An analysis of monte carlointeger programming, simulated annealing, and tabusearch heuristics for solving spatial harvest schedulingproblems.
Forest Science , 292–301 (1999).[32] Sørensen, M. & Stidsen, T. R. Hybridizing inte-ger programming and metaheuristics for solving highschool timetabling. In Proceedings of the 10th interna-tional conference of the practice and theory of automatedtimetabling , 557–560 (2014).[33] Danna, E., Rothberg, E. & Pape, C. L. Exploring relax-ation induced neighborhoods to improve MIP solutions.
Mathematical Programming , 71–90 (2004).[34] Glover, F. Heuristics for integer programming using sur-rogate constraints.
Decision Sciences , 156–166 (1977).[35] Traversa, F. L. & Di Ventra, M. Memcomputing in-teger linear programming. arXiv:1808.09999 (2018).http://arxiv.org/abs/1808.09999v1.[36] Barnhart, C., Johnson, E. L., Nemhauser, G. L., Savels-bergh, M. W. & Vance, P. H. Branch-and-price: Columngeneration for solving huge integer programs. Operationsresearch , 316–329 (1998).[37] Benders, J. F. Partitioning procedures for solving mixed-variables programming problems. Numerische mathe-matik , 238–252 (1962).[38] Hooker, J. N. & Ottosson, G. Logic-based bendersdecomposition. Mathematical Programming , 33–60(2003).[39] Abara, J. Applying integer linear programming to thefleet assignment problem. Interfaces , 20–28 (1989).[40] Kroon, L. et al. The new dutch timetable: The or revo-lution.
Interfaces , 6–17 (2009).[41] Stahlbock, R. & Voß, S. Operations research at containerterminals: a literature update. OR spectrum , 1–52(2008).[42] Melo, M. T., Nickel, S. & Saldanha-Da-Gama, F. Fa-cility location and supply chain management–a review. European journal of operational research , 401–412(2009).[43] Bard, J. F., Binici, C. et al.
Staff scheduling at the unitedstates postal service.
Computers & Operations Research , 745–771 (2003).[44] Traversa, F. L. & Di Ventra, M. Universal memcomput-ing machines. IEEE Trans. Neural Netw. Learn. Syst. , 2702 (2015).[45] Di Ventra, M. & Traversa, F. L. Perspective: Memcom-puting: Leveraging memory and physics to compute effi-ciently. Journal of Applied Physics , 180901 (2018).[46] Traversa, F. L. & Di Ventra, M. Polynomial-time so-lution of prime factorization and np-complete problems with digital memcomputing machines. Chaos: An In-terdisciplinary Journal of Nonlinear Science , 023107(2017).[47] Di Ventra, M. & Traversa, F. L. Absence of chaos indigital memcomputing machines with solutions. Phys.Lett. A , 3255 (2017).[48] Di Ventra, M. & Traversa, F. L. Absence of periodic or-bits in digital memcomputing machines with solutions.
Chaos: An Interdisciplinary Journal of Nonlinear Sci-ence , 101101 (2017).[49] URL .[50] Di Ventra, M. & Traversa, F. L. Self-organizing logicgates and circuits and complex problem solving withself-organizing logic circuits, US patent application No.15/557,641, US patent No. 9,911,080 (2018).[51] Traversa, F. L., Ramella, C., Bonani, F. & Di Ventra,M. Memcomputing NP-complete problems in polynomialtime using polynomial resources and collective states. Science Advances , e1500031 (2015).[52] Traversa, F. L. Collective Computing. In preparation (2018).[53] Di Ventra, M., Traversa, F. L. & Ovchinnikov, I. V. Topo-logical field theory and computing with instantons.
Ann.Phys. (Berlin)
PhysicalReview Applied (2018).[55] Traversa, F. L., Cicotti, P., Sheldon, F. & Di Ventra, M.Evidence of exponential speed-up in the solution of hardoptimization problems. Complexity , 1–13 (2018).[56] Sheldon, F., Cicotti, P., Traversa, F. L. & Di Ventra, M.Stress-testing memcomputing on hard combinatorial op-timization problems.
Preprint arXiv:1807.00107 (2018).http://arxiv.org/abs/1807.00107v1.[57] Sheldon, F., Traversa, F. L. & Di Ventra, M. Taming anon-convex landscape with long-range order.
In prepara-tion .[58] Manukian, H., Traversa, F. L. & Di Ventra, M. Accelerat-ing deep learning with memcomputing.
Neural Networks
In preparation (2019).[61] Sharp, T., , Qian, Z. & Traversa, F. L. Distributing Mem-computing on Graphics Processing Units.
In preparation (2019). APPENDIX function [problem, J, V, VL, xe cg, We, input] = ...Airbus problem generator(xmin cg,xmax cg,n1,n2,n3,N)% airbus 5th problem generatorn = n1+n2+n3; % total number of containersWp = 40000; % Max payloadWe = 120000; % Aircraft weightxe cg = -0.05; % Center of gravity position of the aircraft without payloadSmax 0 = 22000; % Max shear% generate distributions of container massinput = container distribution generator(n1,n2,n3,N);type1 = find(input(:,2)==1);type2 = find(input(:,2)==2);type3 = find(input(:,2)==3);indexM = zeros(n,N);var = 0;for j=1:n1indexM(type1(j),:) = var+(1:N);var = var+N;endfor j=1:n2indexM(type2(j),:) = var+(1:N);var = var+N;endfor j=1:n3indexM(type3(j),1:end-1) = var+(1:N-1);var = var+N-1;end% distance defined in (3) d = zeros(n,N);d(type1,:) = repmat(linspace(-N/2+1/2,N/2-1/2,N)/N,n1,1);d(type2,:) = repmat(linspace(-N/2+1/2,N/2-1/2,N)/N,n2,1);d(type3,1:end-1) = repmat(linspace(-N/2+1,N/2-1,N-1)/N,n3,1);% initialize matrices for IPAineq = [];bineq = [];% constraints on containers Eq. (2a)Aineq = [Aineq;sparse(repmat((1:n1).’,1,N),indexM(type1,:),1,n1,var);sparse(repmat((1:n2).’,1,N),indexM(type2,:),1,n2,var);sparse(repmat((1:n3).’,1,N-1),indexM(type3,1:end-1),1,n3,var)];bineq = [bineq; ones(n1+n2+n3,1)]; % constraints on bins Eq. (2b)Aineq = [Aineq;sparse(1,indexM([type1; type2; type3],1),[ones(n1,1); ones(n2,1)/2; ones(n3,1)],1,...var); sparse(repmat(1:N-2,n1+n2+2*n3,1),[indexM([type1; type2; type3],2:end-1);indexM(type3,1:end-2)],[ones(n1,N-2); ones(n2,N-2)/2; ones(2*n3,N-2)],N-2,var);sparse(1,[indexM(type1,end);indexM(type2,end);indexM(type3,end-1)],[ones(n1,1);ones(n2,1)/2; ones(n3,1)],1,var)];bineq = [bineq; ones(N,1)];% constraint on max weight (2c)Aineq = [Aineq;sparse(1,indexM([type1; type2],:),repmat(input([type1; type2],3),1,N),1,var)+...sparse(1,indexM(type3,1:end-1),repmat(input(type3,3),1,N-1),1,var)];bineq = [bineq; Wp];% constraint on the center of gravity (2d)Aineq = [Aineq;sparse(1,indexM([type1; type2],:),repmat(input([type1; type2],3),1,N).*(d([type1;type2],:)-xmax cg),1,var)+sparse(1,indexM(type3,1:end-1),repmat(input(type3,3),...1,N-1).*(d(type3,1:end-1)-xmax cg),1,var)];bineq = [bineq; -(xe cg-xmax cg)*We];Aineq = [Aineq;sparse(1,indexM([type1; type2],:),repmat(-input([type1; type2],3),1,N).*(d([type1;type2],:)-xmin cg),1,var)+sparse(1,indexM(type3,1:end-1),repmat(-input(type3,3),...1,N-1).*(d(type3,1:end-1)-xmin cg),1,var)];bineq = [bineq; (xe cg-xmin cg)*We];% constraint on the shear curve (2e)for j=1:floor(N/2)Aineq = [Aineq;sparse(1,indexM([type1; type2],1:j),repmat(input([type1; type2],3),1,j),1,var)+...sparse(1,indexM(type3,1:j),[repmat(input(type3,3),1,j-1) input(type3,3)/2],1,var)];bineq = [bineq; Smax 0*j/floor(N/2)];Aineq = [Aineq;sparse(1,indexM([type1; type2],end:-1:end-j+1),repmat(input([type1;type2],3),1,j),1,var)+sparse(1,indexM(type3,end-1:-1:end-j),[repmat(input(type3,...3),1,j-1) input(type3,3)/2],1,var)];bineq = [bineq; Smax 0*j/floor(N/2)];end% objective function (5)f = full(-sparse(1,indexM([type1; type2],:),repmat(input([type1; type2],3),1,N),1,...var)-sparse(1,indexM(type3,1:N-1),repmat(input(type3,3),1,N-1),1,var)); % generate output problem structure can be tested in Matlab using intlinprogproblem.Aeq = [];problem.beq = [];problem.Aineq = Aineq;problem.bineq = bineq;problem.f = f;problem.lb = zeros(var,1);problem.ub = ones(var,1);problem.solver = ’intlinprog’;problem.options = [];problem.intcon = 1:var;% generate indexes to evaluate the center of gravity as% x cg = (sparse(1,J,VL)*X+xe cg*We)./(sparse(1,J,V)*X+We)J = indexM([type1; type2],:);dummy = indexM(type3,1:end-1);J = [J(:); dummy(:)];VL = repmat(input([type1; type2],3),1,N).*d([type1; type2],:);dummy = repmat(input(type3,3),1,N-1).*d(type3,1:end-1);VL = [VL(:); dummy(:)];V = repmat(input([type1; type2],3),1,N);dummy = repmat(input(type3,3),1,N-1);V = [V(:); dummy(:)];function input = container distribution generator(n1,n2,n3,N)w1 = ([(3500-1500)/3*randn(n1*1000,1)+1500; (3500-1500)/3*randn(n1*1000,1)+3500]);w1 = round(w1((w1>1300)&(w1<3700))*20/N);w2 = ([(1800-700)/3*randn(n2*1000,1)+700; (1800-700)/3*randn(n2*1000,1)+1800]);w2 = round(w2((w2>500)&(w2<2000))*20/N);w3 = ([(7000-3200)/3*randn(n3*1000,1)+3200; (7000-3200)/3*randn(n3*1000,1)+7000]);w3 = round(w3((w3>3000)&(w3<7200))*20/N);input = [(1:n1).’ ones(n1,1) w1(randperm(length(w1),n1));n1+(1:n2).’ 2*ones(n2,1) w2(randperm(length(w2),n2));n1+n2+(1:n3).’ 3*ones(n3,1) w3(randperm(length(w3),n3))];% generate output problem structure can be tested in Matlab using intlinprogproblem.Aeq = [];problem.beq = [];problem.Aineq = Aineq;problem.bineq = bineq;problem.f = f;problem.lb = zeros(var,1);problem.ub = ones(var,1);problem.solver = ’intlinprog’;problem.options = [];problem.intcon = 1:var;% generate indexes to evaluate the center of gravity as% x cg = (sparse(1,J,VL)*X+xe cg*We)./(sparse(1,J,V)*X+We)J = indexM([type1; type2],:);dummy = indexM(type3,1:end-1);J = [J(:); dummy(:)];VL = repmat(input([type1; type2],3),1,N).*d([type1; type2],:);dummy = repmat(input(type3,3),1,N-1).*d(type3,1:end-1);VL = [VL(:); dummy(:)];V = repmat(input([type1; type2],3),1,N);dummy = repmat(input(type3,3),1,N-1);V = [V(:); dummy(:)];function input = container distribution generator(n1,n2,n3,N)w1 = ([(3500-1500)/3*randn(n1*1000,1)+1500; (3500-1500)/3*randn(n1*1000,1)+3500]);w1 = round(w1((w1>1300)&(w1<3700))*20/N);w2 = ([(1800-700)/3*randn(n2*1000,1)+700; (1800-700)/3*randn(n2*1000,1)+1800]);w2 = round(w2((w2>500)&(w2<2000))*20/N);w3 = ([(7000-3200)/3*randn(n3*1000,1)+3200; (7000-3200)/3*randn(n3*1000,1)+7000]);w3 = round(w3((w3>3000)&(w3<7200))*20/N);input = [(1:n1).’ ones(n1,1) w1(randperm(length(w1),n1));n1+(1:n2).’ 2*ones(n2,1) w2(randperm(length(w2),n2));n1+n2+(1:n3).’ 3*ones(n3,1) w3(randperm(length(w3),n3))];