Aircraft Loading Optimization: MemComputing the 5^{th} Airbus Problem

Aircraft Loading Optimization: MemComputing the Airbus Problem

Fabio L. Traversa email: MemComputing, Inc., San Diego, CA, 92130 CA
July 26, 2019

On the January 22th 2019, Airbus launched a quantum computing challenge to solve a set of problems relevant for the aircraft life cycle (Airbus challenge web-page). The challenge consists of a set of 5 problems that ranges from design to deployment of aircraft. This work addresses the problem. The formulation exploits an Integer programming framework with a linear objective function and the solution relies on the MemComputing paradigm. It is discussed how to use current MemComputing software MemCPU 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 discussed are possible formulations of the problem utilizing non-linear objective functions, allowing for different optimization schemes implementable in modified MemCPU software, potentially useful for field operation purposes.

I Introduction

Automated computation in the 21st century has reached a paramount role in industry, finance, consumer technology and much more Eiben2015 (); Kephart2003 (). However, the more that computation becomes relevant, the greater the challenges presented to both industry and academia. Solutions to today’s computational limits are provided by more and more sophisticated and high-performing CPUs, GPUs etc. Vestias2014 (), as well as by the insurgence of distributed computing in cloud infrastructures Hashem2015 ().

Nevertheless, the unceasing growth of computing demand is pushing towards completely new architectures as well as new paradigms. This can mean not only “handling more computation,” but also making the computation more efficient. For example, neuromorphic computing based on diverse realizations of artificial neural networks Burr2016 (); Torrejon2017 (); Furber2014 () promises a paradigm shift in machine learning and artificial intelligence.

In this scenario, new computing architectures based on quantum physics and not strictly on algorithmic approaches promise to solve among the most challenging problems ranging from drug discovery to A ILinke2017 (); Bennett2000 (); Ladd2010 (); Raha2007 (). However, it is still not clear if and when reliable quantum computers will show what currently goes under the name of quantum supremacy over classical computing architectures Harrow2017 (); MITtechrev (); Dyakonov2019 (); Tang2018 ().

However, companies and academia are already looking at applications for future quantum computers Accenture_report (); MITtechrev (); AribusURL (). I focus here on the challenge recently launched by Airbus AribusURL () consisting of 5 problems related to the life cycle of an aircraft, i.e., ranging from the design to the deployment of an aircraft. The general requirements for the challenge are that for each problem a suitable algorithm implementable on a quantum computer should be developed, tested on a simulated quantum computer, and finally assessed in terms of the performance at scale.

I consider in this work the 5th Airbus problem AribusURL (): the Aircraft Loading Optimization (ALO) problem. The ALO requires optimizing the placement of containers of different sizes and weights in an aircraft subject to limitations on maximum weight allowed, maximum tolerable shear and center of gravity. However, instead of developing on algorithm for future quantum computers, I present a formalization of the problem using Integer Programming (IP) framework Schrijver1998 (). IP is a problem formulation widely used in industry and academia consisting of an objective function to be minimized over a set of variables. Variables may be constrained to take binary, integer or continuous values. Moreover, the objective function may be subject to further constraints in the form of linear inequalities among variables. Finally, depending on the nature of the objective function, we can have different IP problems such as linear, quadratic or other non-linear type.

The formulation I consider here covers all requirement of the ALO statement and leads to an IP problem involving binary variables only and having two objective functions to be minimized, one linear and the other non-linear. I therefore propose two different solution strategies to solve this multi-objective IP problem.

The binary IP (BP) problem is NP-complete, sometimes also known as 0-1 linear programming, and is one of Karp’s 21 NP-complete problems complexity_bible (). Therefore, it is not surprising that Airbus poses a challenge about a problem whose natural formulation leads to a special case of on NP-complete problem. Several general-purpose open source GLPK (); Forrest2018 (); Gleixner2018 (); Ralphs2017 (); Berkelaar2004 () and commercial solvers Gurobi (); CPLEX (); FICO (); MATLAB () have been developed to solve IP problems. However, these can fail when the problems are particularly hard and therefore specialized solvers optimized to exploit specific structures for ILP have also been developed applegate2006concorde (); Floudas2005 (); Boston99 (); sorensen2014hybridizing ().

Figure 1: Sketch of the Airbus Aircraft Loading Optimization problem. The aircraft is to be loaded by selecting a portion of the available payload. The containers can have different shapes. For this ALO, 3 different shapes are considered numbered from 1 to 3 as in the figure. The maximum weight allowed is limited. The center of gravity of the system aircraft + containers must remain within the limit and the shear curve must also be limited to remain under a give shear limit curve .

IP problems can be approached employing different techniques. A classification of those includes heuristics Danna2004 (); Boston99 (); Glover1977 (); Traversa2018a () and exhaustive algorithms Schrijver1998 (); barnhart1998branch (); benders1962partitioning (); hooker2003logic () also called “complete algorithms.” The complete algorithm for ILP that is most commonly employed is a combination of cutting planes and branch-and-bound, also known as the branch-and-cut algorithm Schrijver1998 (). All these solvers have demonstrated varied degrees of success on a variety of IP problems abara1989applying (); kroon2009new (); stahlbock2008operations (); melo2009facility (); bard2003staff (), however scaling properties for BP problems such as the one proposed by Airbus can easily show exponential blowing-up when applying one of those methods complexity_bible ().

In this work I discuss the general purpose approach to the solution of the Aibus ALO problem based on the memcomputing paradigm introduced in UMM (). Memcomputing (computing with and in memory Di_Ventra2018 ()) approach is neither based on stochastic search nor on trial and error strategy. In addition, memcomputing does not use a set of instructions recursively employed to find solutions to the problem at hand. Therefore, we can classify the memcomputing approach as non-algorithmic Di_Ventra2018 ().

The memcomputing approach relies on embedding a given problem into an electronic circuit, which represents a possible realization of a memcomputing machine (MM) UMM (); DMM2 (); Di_Ventra2018 (). These electronic circuits, when designed in order to satisfy a few properties DMM2 (), naturally relax to an equilibrium that represents the solution of the problems at hand DMM2 (); Di_Ventra2018 (); no-chaos (); noperiod ().

In order to approach the Airbus ALO problem I use the MM realization described in Traversa2018a () employing self-organizing algebraic gates (SOAGs). The SOAG-based MM can be efficiently simulated in software and in this work I use the software developed at MemComputing Inc. MemWeb () named MemCPU.

In section II I briefly describe the Airbus ALO problem. In section III the BP formulation of the ALO problems is discussed. In section IV a quick overview of memcomputing and MemCPU is provided, while in section V are discussed both numerical results and scaling properties of MemCPU applied to the ALO problem.

Ii Airbus Aircraft Loading Optimization Problem

The ALO problem proposed by Airbus is sketched in Fig. 1. Following the problem statement we need to optimize the placement of a subset of containers picked from an available payload formed by different containers. Each container is described by the triplet with the identification number, the mass and the size of the container.

The aircraft has available positions for standard cargo as depicted in Fig. 1. The size can be of 3 different types: size 1 occupies one position; size 2 occupies half a position and may share the position with another size 2 container; size 3 occupies 2 positions.

The problem requires maximizing the mass of the carried freight for the flight such that:

  1. the containers must be placed consistently with size and positions available

  2. the mass of the carried freight does not exceed the maximum payload capacity of the aircraft

  3. the center of gravity position of the carried freight + aircraft must be within the limits

  4. the shear curve defined by the container distribution must be bounded by a given maximal shear curve

  5. while maximizing the mass, also the position of the center of gravity of the carried freight + aircraft should be optimized to be as close as possible to a target center of gravity .

These constraints are sketched also in Fig. 1. While the constraints (A) and (B) do not need additional details or explanations, the others do. In (C) the center of gravity is defined by the the equilibrium of all forces (container + aircraft) in the gravity axis direction (see force in Fig. 1). Each container is supposed to have uniform mass, therefore the weight force can be considered as the weight of a point located at the center of the container. Finally, the mass and location of the center of gravity of the empty aircraft are provided.

The constraint (D) requires a consistent definition of the shear curve. The shear curve is defined by


where is the length of the loading region of the aircraft and is the position relative to the center of the loading zone (see Fig. 1). In the statement of the problem, the maximum shear curve can be either linear and symmetric with respect to , or asymmetric or some non-linear function of .

Finally, (E) represents an extra optimization requirement and not really a constraint. In fact, from the Airbus statement, what should be primarily optimized is the mass of the carried freight, then optimize the distribution such that the resulting center of gravity is as close as possible to a target.

Iii Integer Programming Formulation

A natural way to formulate mathematically the problem posed in section II is using Integer Programming framework Schrijver1998 ().

Let us start assigning a binary variable to each container located in the position of the aircraft. Since each container can be located in at most one position in the aircraft we constrain these variables requiring


On the other hand, we have the constraint on the size and the number of available spots on the aircraft that can be also expressed as a non-linear inequality. In order to formalize this constraint, let us further characterize the variable . We consider here, for simplicity, only the 3 sizes defined in the Airbus statement. However, this formulation can be easily extended to more complicated sizes and bin distributions in the aircraft. Let us define and as three sets of indexes such that . Moreover, if then . Therefore, and regroup the variables in sets of variables corresponding to the same container size.

Using this index characterization, the constraint (A) on the sizes and bins can be expressed as


It is easy to prove that these inequalities guarantee that for each bin there is no overlapping of containers except for containers of size 2 for which two of them can occupy the same bin (the coefficient allows the overlapping). The case of the size 3 (containers that occupy two consecutive bins) is enforced by the term in which two consecutive containers of size 3 appear in the same constraint; therefore a container of size 3 occupies 2 bins if selected. It is also worth noticing that containers of size 3 have that ranges only from 1 to because, since they occupy two bins, there ore only possible locations for them.

The constraint (B) can be trivially described by


The constraint (C) on the center of gravity requires the definition of signed distance from for a container of size located in the bin . Considering equally distributed bins around (this assumption is the same reported in the Airbus statement but can be easily relaxed), we have that


Using this distance definition, the center of gravity of the loaded aircraft can be evaluated as


where and are respectively the mass and the center of gravity of the empty aircraft.

Eq. (4) allows the formalization of the constraint (C) through the inequalities


that express and respectively.

Regarding the constraint (D), the shear curve can be easily calculated at each bin location and set smaller than :


which expresses the shear curve defined in Eq. (II) at the centers of the bins taking into account the different sizes of the containers.

Finally we can define the objective function of the ALO problem as


whose minimization over subject to the collection of constraints (III) provides the distribution and the maximum mass of the carried freight satisfying the constraints (A)-(D) of section II.

iii.1 Center of Gravity Optimization

We briefly discuss in this section the further optimization required in (E) of section II. This can be handled with the following scheme. Once the problem (5) is solved, we obtain the maximum mass of the carried freight as


where is the solution of the minimization of (5) subject to (III). We therefore define a new IP problem where we include the constraints (IIIa,IIIb,IIIc,IIIe) but exclude for now (IIId). Instead we include the extra constraint


where is a tolerance parameter that can be used in case we are interested in a slightly lighter carried freight but better center of gravity.

We then define the new objective function


with defined by Eq. (4).

The target now is to find the assignment of that minimizes Eq. (8) under the constraints (IIIa,IIIb,IIIc,IIIe) and (7). This can be achieved either using a sequence of linear IP problems in which we reintroduce the constraint (IIId) and we iteratively shrink the around , or by implementing it directly as a non-linear IP problem with objective function Eq. (8).

Iv MemComputing Approach

The memcomputing approach to IP problems is based on the concept of Self-Organizing Algebraic Gates (SOAGs) introduced in Traversa2018a (). SOAG is a novel circuit design developed at MemComputing, Inc. MemWeb () by the author of this work and is part of a class of self-organizing circuits as Self-Organizing Logic Gates (SOLGs) introduced in DMM2 (); patentSOLC ().

Both SOLGs and SOAGs are building blocks for practical realizations of Universal Memcomputing Machines (UMM) Di_Ventra2018 (); UMM (); traversaNP (), with digital input-output (Digital MM, DMM) DMM2 ().

Figure 2: Reprint from Traversa2018a (). Sketch of a Self-Organizing Algebraic Gate. All terminals allow a superposition of incoming and outgoing signals from the surrounding circuit. The central unit processes the signals in order to satisfy a linear algebraic relation consistent with the requirement of the “out” terminal. The self-organization is enforced by the Dynamic Correction Modules that read voltages from all terminals and inject a current to the appropriate terminal as long as the algebraic relation is not satisfied.

Figure 3: Reprint from Traversa2018a (). Sketch of a Self-Organizing Algebraic Circuit (SOAC). SOAGs are connected together in an architecture that directly maps the IP into the SOAC.

The SOAG is designed to self-organize toward an algebraic relation representing a linear inequality among variables. In this work, the SOAGs are designed to satisfy linear relations between binary variables (see Fig. 2) since the ALO problem formulation involves binary variables only.

By connecting together SOAGs we form a Self-Organizing Algebraic Circuit (SOAC), see Fig. 3. The SOAC collectively self-organizes in order to satisfy the relations embedded in the gates. In this way, it is trivial to embed an IP problem directly into the SOAC. Each inequality in (III) or (7) can be mapped directly into a SOAG. The objective functions (5) can be easily reformulated as an extra linear inequality


where represents a threshold for the maximum carried freight that is dynamically changed each time a feasible solution is found in order to find solutions closer and closer to the global minimum of the problem. On the other hand, the objective function (8) can be implemented as a pair of linear inequalities of the form


where is an extra parameter that again can be dynamically changed to find solutions of increasing quality, each time closer to the global optimum. and are linear functions of defined by substituting (4) in (8). This also leads to defining as the threshold of the distance of the actual center of gravity from the target .

The problem formulated and embedded in a SOAC as described can be efficiently handled by either actually building the electronic circuit or just by simulating it since it involves only standard (non-quantum) electronic components. Simulating means solving differential equations of the form


with appropriate initial conditions and where and are non-linear functions of the voltages (), currents () and extra internal state variables () characterizing the circuit. In this picture a configuration of the voltages () at a given time represents an actual assignment to the variables of the IP problem by reading the voltages though thresholds: if a voltage at a node of the circuit is above the threshold, then it corresponds to a logical 1, and otherwise corresponds to a logical 0. On the other hand, the transition function of these machines (namely the function that maps input to output) is physical (analog) and takes full advantage of the collective state of the system to process information collective (); traversaNP (); DMM2 (). However, despite its analog nature, the mapping of voltages into binary variable through thresholds allows efficient size scaling for these machines, avoiding the bottleneck related to the precision of writing and reading inputs and outputs. Finally, it is worth noticing here that the memcomputing approach used to solve binary as presented in Traversa2018a () for the case of IP problems, does not provide proof of optimality for a given solution, nor does it detect the infeasibility of a problem.

In this work we consider the simulation of the SOAC implemented in the software MemCPU developed at MemComputing, Inc.MemWeb () and available also as Software at a Service.

The working principle of SOAGs and SOACs, i.e., self-organization of voltages and currents of an electronic circuit in order to satisfy algebraic relations, is enforced by both active and passive electronic elements with and without memory Traversa2018a (); DMM2 (); Di_Ventra2018 (). The features of DMMs have been investigated, and interesting properties emerge from a correct design such as long-range order correlations and topological robustness topo (); Bearden2018 (). If mathematical requirements described in DMM2 () are fulfilled, then persistent oscillations or chaos are avoided no-chaos (); noperiod (). Self-organizing digital circuits (i.e. both SOLCs and SOACs), being a realization of DMMs can in principle solve efficiently (i.e. with polynomial resources) complex combinatorial optimization problems, given that mathematical features on the constitutive equation (12) are satisfied UMM (); DMM2 (); Di_Ventra2018 (). Self-organizing digital circuits (i.e. both SOLCs and SOACs) have also been proved to efficiently handle a variety of combinatorial optimization problems ranging from maximum satisfiability (MAXSAT) Traversa2018 (); Sheldon2018 () to quadratic unconstrained binary optimization (QUBO) spinglass () and from IP Traversa2018a () to training of neural networks AcceleratingDL ().

V Scaling Results

In order to assess the scaling properties of solving the Airbus ALO problem employing MemCPU we need to generate a set of meaningful benchmarks at different values of and . Since we have no information about the actual values of and and not even on the ratio between number of containers at different sizes or typical mass distribution of the containers available, we use the sample data set provided by Airbus (see Table 1 for the data set, or it can also be found from AribusURL ()) to extrapolate. Notice that the data set provided from Airbus is just “for illustration and for testing the algorithm” AribusURL ().

, linear, symmetric
1 1 2134
2 1 3455
3 1 1866
4 1 1699
5 1 3500
6 1 3332
7 1 2578
8 1 2315
9 1 1888
10 1 1786
11 1 3277
12 1 2987
13 1 2534
14 1 2111
15 1 2607
16 1 1566
17 1 1765
18 1 1946
19 1 1732
20 1 1641
21 2 1800
22 2 986
23 2 873
24 2 1764
25 2 1239
26 2 1487
27 2 769
28 2 836
29 2 659
30 2 765
Table 1: Airbus ALO data set AribusURL ().
Figure 4: Histograms of container mass from Airbus data set reported in Table 1.

In order to generate benchmarks for different sizes we need to generate a distribution of containers. One of the possible ways is to try to produce a distribution similar to that from the sample data set. In Fig. 4 the distributions of container masses from the data set in Table 1 is reported. Even though there are not enough data to certainly determine the distribution shape, it is reasonable to assume that both distributions can be recreated from a bimodal Gaussian distribution, cut off at certain boundaries.

Figure 5: Monte Carlo distribution of the objective (5) versus MemCPU parameters MemWeb (). The ALO problem is defined with the parameters in Table 1 and . Number of Monte Carlo iterations = 1200, simulation time = 8 and number of Markov Chains = 8.

In Fig. 11 I have included Matlab code that generates ALO problems for any and (it is also available a converter from Matlab to .mps format at LuongMatlab ()). The distribution of containers is generated in the subfunction container_distribution_generator. A bimodal Gaussian distribution of containers is generated for each size type. The ranges are chosen in such a way that the distributions in Fig. 4 can be qualitatively reproduced. It is also worth noticing that the masses are scaled depending on . In fact, maintaining the same weight of the airplane and the same limit on the maximum weight, but increasing the number of bins, requires for consistency that each bin be scaled in size, and therefore we scale all container masses accordingly.

v.1 MemCPU results

Time Constants Static Constants
Time Const. 1 0.79293 Static Const. 1 0.86027
Time Const. 2 0.28953 Static Const. 2 0.43655
Time Const. 3 0.19784 Static Const. 3 0.24784
Time Const. 4 0.85076 Static Const. 4 0.68352
Time Const. 5 0.42802 Static Const. 5 0.37106
Time Const. 6 0.39473 Static Const. 6 0.44440
Time Const. 7 0.76838 Static Const. 7 0.26562
Static Const. 8 0.21645
Static Const. 9 1.00000
Static Const. 10 0.88580
Table 2: Parameters for MemCPU MemWeb () extracted from the distribution of Fig. 5.
Figure 6: Configuration of forces to evaluate how close the center of gravity can approach the target (see Table 1), without restriction on the minimum selected payload () and with maximum selected payload (). For the maximum selected payload, the maximum density of (container mass per unit of aircraft length) is assumed, consistent with the data set in Table 1 and the subfunction container_distribution_generator in Fig. 11.

As discussed in section IV, MemCPU is an emulator of an electronic circuit with a given IP problem embedded within. However, in order to be efficient in solving the problem at hand, MemCPU needs a proper set of parameters characterizing electronic elements such as resistors, capacitors, etc. These parameters do not depend on the size of the problem (they are scale free) but depend on the structure of the problem. MemComputing currently provides a Monte Carlo routine that evaluates the distribution of the objective function of the IP problem as a function of the parameters tuning (); however future releases will provide a predictor routine for parameters that will avoid the Monte Carlo-based tuning tuning (). The objective function distribution is reported in Fig. 5. This distribution has been computed for the problem generated using the code in Fig. 11 using as input the container distribution in Table 1 and . The choice of is only for tuning purposes because this restricts the possible feasible solutions of the problem making it “harder.” In fact, it is easy to estimate the limit of the center of gravity for any selection of containers () from the available payload. In Fig. 6, the forces in the center of gravity evaluation (Eq. 4) are summarized for the limit of containers of all equal weight corresponding to the maximum weight allowed by the maximum shear curve. The value is and I chose that largely reduces the number of feasible solutions. This choice for leads to minimal feasible solutions that have objective far from maximum allowed carried freight , however, used for tuning is a good choice because largely reduces the number of feasible solutions. This choice helps to produce a sharper distribution of objective function values versus parameters which also is useful for solving the problem in section III.1.

From the objective distribution of Fig.5, a set of the parameters for MemCPU can be easily extracted and used for running any other instance of the ALO problem. However, extracting a good set of parameters 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 returning the best objective is selected. If there is more than one set of parameters with the same objective for the same Markov chain, the choice is made based upon the one that has the best average of the objectives related to the closest 10 sets of parameters within the same Markov chain. This super-selection rule comes from the fact that each MemCPU run starts with random initial conditions for voltages and currents, and therefore there could be a fluctuation in the output objective due to the finite (and short) simulation time tuning (). Once we have the best parameter choice per Markov chain, we can run the same problem multiple times with random initial conditions for the same set of parameters and select the set of parameters that statistically arrives earlier to the best objective in a given time out. In this way, I extracted the set of parameters reported in Table 2.

Figure 7: MemCPU run time versus to find a configuration of containers with total mass larger that than 99.9% of the maximum mass allowed on the aircraft or of the total payload if smaller than the maximum mass allowed. Dots connected by solid lines are the average time for 100 different ALO instances generated using the code in Fig. 11. Dashed curves are the scaling relation (13). Runs have been carried out on an Intel(R) Xeon(R) Gold 6138 CPU @ 2.00GHz.
Figure 8: MemCPU run time versus to find a configuration of containers with total mass larger that than 99.9% of the maximum mass allowed on the aircraft or of the total payload if smaller than the maximum mass allowed. Dots connected by solid lines are the average time for 100 different ALO instances generated using the code in Fig. 11. Dashed curves are the scaling relation (13). Runs have been carried out on an Intel(R) Xeon(R) Gold 6138 CPU @ 2.00GHz.

Using the parameters in Table 2, MemCPU has been tested on a set of ALO instances generated for different and using the code in Fig. 11. The created benchmark consists of 100 ALO instances for each pair . Each instance has a different container mass distribution randomly generated by the subfunction container_distribution_generator of Fig. 11. For each instance, the container sizes were , and with appropriate rounding, in order to have integer and .

Since the value of the objective related to the global minimum of the ALO problem is not known and MemCPU does not provide proof of optimality, I have set a very tight threshold to assess the quality of the solution: accept the solution if the objective function is smaller than if or smaller than otherwise, with . This is equivalent to requiring a container selection that is larger than 99.9% of either the maximum weight allowed on the aircraft or of the total payload if the latter is smaller than the former. In Fig. 7 and 8 the time to find the solution for is reported as a function of (Fig. 7) or as a function of number of nonzero entries in the matrix defined by constraints (III) (Fig. 8) for different ratios .

The results as a function of are useful for properly assessing the MemCPU scaling. Indeed, MemCPU associates to each constraint a SOAG with a number of terminals equal to the terms in the constraint. For each terminal there is a dynamic correcting module (DCM) to simulate Traversa2018a (); DMM2 (). Therefore, the number of DCMs that MemCPU simulates is equal to , and it is natural to measure the complexity in terms of that is at the same time also a measure of the size of an IP problem. The interpolation curve on - scale in the plane -time provides the following relation:


where the dependence on is valid for . The equation (13) shows a sub-quadratic scaling of the time to solution versus . From the (III) it is easy to realize that and, substituting this relation in equation (13), the scaling as a function of and can be recovered.

It is worth discussing the dependence of Eq. (13) on . The first thing that can be noticed is that the dependence of on is very weak. The coefficient 0.11 in (13) makes the exponent range from 1.26 to 1.58 for . On top of the weak dependence, there is compensation from the stronger and negative dependence of the prefactor on . However, increasing , an evident dependence of the exponent on disappears, and it saturates for . This is not surprising behavior because for larger the problem becomes “easier” since there are many more choices of container selections that respect the constraints, and the burden in the calculation depends only on the size of the ALO. We do not show here numerical experiments to support the latter claim because considering much larger the seems meaningless in practical cases.

Figure 9: MemCPU run time versus to find a configuration of containers that a) has total mass larger that than 99.8% of the maximum mass allowed on the aircraft or of the total payload if smaller than the maximum mass allowed and b) optimizes the center of gravity of the aircraft. Dots connected by solid lines are the average time for 100 different ALO instances generated using the code in Fig. 11. Dashed curves are the scaling relation (13). Runs have been carried out on an Intel(R) Xeon(R) Gold 6138 CPU @ 2.00GHz.

Let us consider now the ALO problem described in (section III.1). In order to handle this problem just using a non-modified version of MemCPU, I have created a sequence of linear IP problems that converges to the solution of the non linear problem of (section III.1). This sequence of problems can be easily created making use of the code of Fig. 11 as follows:

  1. Set and for some initial values, for example, the ones given in Table 1.

  2. Generate the problem for a given and .

  3. Add to the matrix of the constraints the objective function of the problem as an extra constraint. As the right hand side of this constraint set where if or otherwise.

  4. Substitute the objective function with a null objective function and solve the problem.

  5. If a solution to the problem is found, evaluate by means of J, VL, V, xe_cg and We as reported in the code of figure 11

  6. if set otherwise set for some . Go back to b.

Figure 10: Center of gravity of the system aircraft + selected payload with total carried mass larger than of . In the inset, the mean of the total mass generated by the subfunction container_distribution_generator of Fig. 11 and versus is plotted.

In Fig. 9 the time to solution is reported for selecting the payload that maximizes the mass ( has been used in this case) and at the same time provides the center of gravity as close as possible to the target, for different value of and . The dashed curves that interpolate this version of the ALO results as a shift of the equation (13) where the coefficient is replaced by . Finding the same scaling properties as in the other problem is not surprising, since solving this ALO problem follows the steps a.-f., and so is a cascade of problems similar to the ones solved in Fig. 7.

Finally, it is worth noticing where the optimized centers of gravity are located. In Fig. 10, we can see that for the centers of gravity approach the first limit of Fig. 6. This is not surprising since in this case, because of the distribution of containers, (see inset of Fig. 10). On the other hand, if , (see inset of Fig. 10) and the optimized centers of gravity approach the second limit of Fig. 6.

v.2 Further improvements

We briefly discuss here a few improvements that can be made in order to have an even more efficient time-to-solution for the Airbus ALO problem computed through MemCPU. However, since they do not qualitatively affect the scaling assessed in this work, the implementation of these is beyond the current scope.

The first improvement concerns the initial values of the objective function (5). MemCPU searches solutions that have an objective equal to or lower than a target. Once found it decreases the objective and searches for the next solution (see section IV and Traversa2018a ()). In this work I have used infinite as the initial objective. With this choice, MemCPU provides many solutions at lower and lower objectives, progressing from the larger to smaller values. Even though this convergence is usually fast, this can be improved even more by giving as a starting value for the objective function .

Another improvement that scales the time to solution down by a few orders of magnitude is running MemCPU on GPUs rather than CPUs SharpGPU (). In fact, the computation performed by MemCPU is nothing else than a circuit simulation that can efficiently distributed on GPUs SharpGPU ().

Finally, the search for maximal carried freight that at the same time optimizes the center of gravity can be accelerated by implementing directly the nonlinear objective function (8), avoiding the sequence of ALO problems discussed in section V.1 and solving only one ALO problem.

Vi Conclusions

In summary, I have shown how to employ MemComputing through the software MemCPU to tackle the 5th Airbus problem efficiently without a quantum computer.

In order to use MemCPU for the 5th Airbus problem, which is an aircraft loading optimization problem, an IP problem has been formulated that includes all constraints required by Airbus with no exception and no approximation.

The scaling properties assessed in this work show sub-quadratic scaling of MemCPU as a function of the size of the problem measured by the number of nonzero elements of the constraint matrix of the associated IP problem.

The scaling properties of MemCPU allow efficient solutions of large ALO problems and it represents a solution ready to be deployed on the field.

Finally, I have also discussed extra improvements that can be made in order to further (and substantially) accelerate the time to solution for the ALO problem.

Vii Acknowledgments

I thank the MemComputing team for useful discussions and in particular Dr. Tristan Sharp for the revision of the manuscript. This work has been supported by MemComputing, Inc.

Viii References


  • (1) Eiben, A. & Smith, J. Introduction to Evolutionary Computing (Springer Berlin Heidelberg, 2015).
  • (2) Kephart, J. & Chess, D. The vision of autonomic computing. Computer 36, 41–50 (2003).
  • (3) Vestias, M. & Neto, H. Trends of CPU, GPU and FPGA for high-performance computing. In 2014 24th International Conference on Field Programmable Logic and Applications (FPL) (IEEE, 2014).
  • (4) Hashem, I. A. T. et al. The rise of “big data” on cloud computing: Review and open research issues. Information Systems 47, 98–115 (2015).
  • (5) Burr, G. W. et al. Neuromorphic computing using non-volatile memory. Advances in Physics: X 2, 89–124 (2016).
  • (6) Torrejon, J. et al. Neuromorphic computing with nanoscale spintronic oscillators. Nature 547, 428–431 (2017).
  • (7) Furber, S. B., Galluppi, F., Temple, S. & Plana, L. A. The SpiNNaker project. Proceedings of the IEEE 102, 652–665 (2014).
  • (8) Linke, N. M. et al. Experimental comparison of two quantum computing architectures. Proceedings of the National Academy of Sciences 114, 3305–3310 (2017).
  • (9) Bennett, C. H. & DiVincenzo, D. P. Quantum information and computation. Nature 404, 247–255 (2000).
  • (10) Ladd, T. D. et al. Quantum computers. Nature 464, 45–53 (2010).
  • (11) Raha, K. et al. The role of quantum mechanics in structure-based drug design. Drug Discovery Today 12, 725–731 (2007).
  • (12) Harrow, A. W. & Montanaro, A. Quantum computational supremacy. Nature 549, 203–209 (2017).
  • (13) MIT Technology Review.
  • (14) Dyakonov, M. When will useful quantum computers be constructed? not in the foreseeable future, this physicist argues. here’s why: The case against: Quantum computing. IEEE Spectrum 56, 24–29 (2019).
  • (15) Tang, E. A quantum-inspired classical algorithm for recommendation systems. arXiv:1807.04271 (2018). eprint
  • (16) Accenture report.
  • (17) Airbus challenge press-release; Airbus challenge web-page.
  • (18) Schrijver. Theory of Linear Integer Programming (John Wiley & Sons, 1998).
  • (19) Garey, M. R. & Johnson, D. S. Computers and Intractability; 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 Optimization Suite 6.0. ZIB-Report 18-26, Zuse Institute Berlin (2018). URL
  • (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
  • (25) URL
  • (26) URL
  • (27) URL
  • (28) URL
  • (29) Applegate, D., Bixby, R., Chvatal, V. & Cook, W. Concorde tsp solver (2006).
  • (30) Floudas, C. A. & Lin, X. Mixed integer linear programming in process scheduling: Modeling, algorithms, and applications. Annals of Operations Research 139, 131–162 (2005).
  • (31) Boston, K. & Bettinger, P. An analysis of monte carlo integer programming, simulated annealing, and tabu search heuristics for solving spatial harvest scheduling problems. Forest Science 45, 292–301 (1999).
  • (32) Sørensen, M. & Stidsen, T. R. Hybridizing integer programming and metaheuristics for solving high school timetabling. In Proceedings of the 10th international conference of the practice and theory of automated timetabling, 557–560 (2014).
  • (33) Danna, E., Rothberg, E. & Pape, C. L. Exploring relaxation induced neighborhoods to improve MIP solutions. Mathematical Programming 102, 71–90 (2004).
  • (34) Glover, F. Heuristics for integer programming using surrogate constraints. Decision Sciences 8, 156–166 (1977).
  • (35) Traversa, F. L. & Di Ventra, M. Memcomputing integer linear programming. arXiv:1808.09999 (2018). eprint
  • (36) Barnhart, C., Johnson, E. L., Nemhauser, G. L., Savelsbergh, M. W. & Vance, P. H. Branch-and-price: Column generation for solving huge integer programs. Operations research 46, 316–329 (1998).
  • (37) Benders, J. F. Partitioning procedures for solving mixed-variables programming problems. Numerische mathematik 4, 238–252 (1962).
  • (38) Hooker, J. N. & Ottosson, G. Logic-based benders decomposition. Mathematical Programming 96, 33–60 (2003).
  • (39) Abara, J. Applying integer linear programming to the fleet assignment problem. Interfaces 19, 20–28 (1989).
  • (40) Kroon, L. et al. The new dutch timetable: The or revolution. Interfaces 39, 6–17 (2009).
  • (41) Stahlbock, R. & Voß, S. Operations research at container terminals: a literature update. OR spectrum 30, 1–52 (2008).
  • (42) Melo, M. T., Nickel, S. & Saldanha-Da-Gama, F. Facility location and supply chain management–a review. European journal of operational research 196, 401–412 (2009).
  • (43) Bard, J. F., Binici, C. et al. Staff scheduling at the united states postal service. Computers & Operations Research 30, 745–771 (2003).
  • (44) Traversa, F. L. & Di Ventra, M. Universal memcomputing machines. IEEE Trans. Neural Netw. Learn. Syst. 26, 2702 (2015).
  • (45) Di Ventra, M. & Traversa, F. L. Perspective: Memcomputing: Leveraging memory and physics to compute efficiently. Journal of Applied Physics 123, 180901 (2018).
  • (46) Traversa, F. L. & Di Ventra, M. Polynomial-time solution of prime factorization and np-complete problems with digital memcomputing machines. Chaos: An Interdisciplinary Journal of Nonlinear Science 27, 023107 (2017).
  • (47) Di Ventra, M. & Traversa, F. L. Absence of chaos in digital memcomputing machines with solutions. Phys. Lett. A 381, 3255 (2017).
  • (48) Di Ventra, M. & Traversa, F. L. Absence of periodic orbits in digital memcomputing machines with solutions. Chaos: An Interdisciplinary Journal of Nonlinear Science 27, 101101 (2017).
  • (49) URL
  • (50) Di Ventra, M. & Traversa, F. L. Self-organizing logic gates and circuits and complex problem solving with self-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 polynomial time using polynomial resources and collective states. Science Advances 1, e1500031 (2015).
  • (52) Traversa, F. L. Collective Computing. In preparation (2018).
  • (53) Di Ventra, M., Traversa, F. L. & Ovchinnikov, I. V. Topological field theory and computing with instantons. Ann. Phys. (Berlin) 1700123 (2017).
  • (54) Bearden, S. R., Manukian, H., Traversa, F. L. & Di Ventra, M. Instantons in self-organizing logic gates. Physical Review Applied 9 (2018).
  • (55) Traversa, F. L., Cicotti, P., Sheldon, F. & Di Ventra, M. Evidence of exponential speed-up in the solution of hard optimization problems. Complexity 2018, 1–13 (2018).
  • (56) Sheldon, F., Cicotti, P., Traversa, F. L. & Di Ventra, M. Stress-testing memcomputing on hard combinatorial optimization problems. Preprint arXiv:1807.00107 (2018). eprint
  • (57) Sheldon, F., Traversa, F. L. & Di Ventra, M. Taming a non-convex landscape with long-range order. In preparation .
  • (58) Manukian, H., Traversa, F. L. & Di Ventra, M. Accelerating deep learning with memcomputing. Neural Networks 110, 1–7 (2019).
  • (59) /19618-mps-format-exporting-tool.
  • (60) Pederson, E., Foertsch, J., Qian, Z. & Traversa, F. L. Tuning and prediction of optimal parameters for algorithm configuration. In preparation (2019).
  • (61) Sharp, T., , Qian, Z. & Traversa, F. L. Distributing Memcomputing on Graphics Processing Units. In preparation (2019).


function [problem, J, V, VL, xe_cg, We, input] = …

x                                    Airbus_problem_generator(xmin_cg,xmax_cg,n1,n2,n3,N)

% airbus 5th problem generator

n = n1+n2+n3; % total number of containers

Wp = 40000; % Max payload

We = 120000; % Aircraft weight

xe_cg = -0.05; % Center of gravity position of the aircraft without payload

Smax_0 = 22000; % Max shear

% generate distributions of container mass

input = 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:n1

x    indexM(type1(j),:) = var+(1:N);

x    var = var+N;


for j=1:n2

x    indexM(type2(j),:) = var+(1:N);

x    var = var+N;


for j=1:n3

x    indexM(type3(j),1:end-1) = var+(1:N-1);

x    var = var+N-1;


% distance defined in (III)

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 IP

Aineq = [];

bineq = [];

% constraints on containers Eq. (IIIa)

Aineq = [Aineq;




bineq = [bineq; ones(n1+n2+n3,1)];

% constraints on bins Eq. (IIIb) Aineq = [Aineq; sparse(1,indexM([type1; type2; type3],1),[ones(n1,1); ones(n2,1)/2; ones(n3,1)],1,… x    var); sparse(repmat(1:N-2,n1+n2+2*n3,1),[indexM([type1; type2; type3],2:end-1); x    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); x    ones(n2,1)/2; ones(n3,1)],1,var)]; bineq = [bineq; ones(N,1)]; % constraint on max weight (IIIc) Aineq = [Aineq; sparse(1,indexM([type1; type2],:),repmat(input([type1; type2],3),1,N),1,var)+… x    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 (IIId) Aineq = [Aineq; sparse(1,indexM([type1; type2],:),repmat(input([type1; type2],3),1,N).*(d([type1; x    type2],:)-xmax_cg),1,var)+sparse(1,indexM(type3,1:end-1),repmat(input(type3,3),… x    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; x    type2],:)-xmin_cg),1,var)+sparse(1,indexM(type3,1:end-1),repmat(-input(type3,3),… x    1,N-1).*(d(type3,1:end-1)-xmin_cg),1,var)]; bineq = [bineq; (xe_cg-xmin_cg)*We]; % constraint on the shear curve (IIIe) for j=1:floor(N/2) x    Aineq = [Aineq; x    sparse(1,indexM([type1; type2],1:j),repmat(input([type1; type2],3),1,j),1,var)+… x    sparse(1,indexM(type3,1:j),[repmat(input(type3,3),1,j-1) input(type3,3)/2],1,var)]; x    bineq = [bineq; Smax_0*j/floor(N/2)]; x    Aineq = [Aineq; x    sparse(1,indexM([type1; type2],end:-1:end-j+1),repmat(input([type1; x        type2],3),1,j),1,var)+sparse(1,indexM(type3,end-1:-1:end-j),[repmat(input(type3,… x        3),1,j-1) input(type3,3)/2],1,var)]; x    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,… x    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 intlinprog problem.Aeq = []; problem.beq = []; problem.Aineq = Aineq; problem.bineq = bineq; problem.f = f; = 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)); x    n1+(1:n2).’ 2*ones(n2,1) w2(randperm(length(w2),n2)); x    n1+n2+(1:n3).’ 3*ones(n3,1) w3(randperm(length(w3),n3))];
Figure 11: Matlab code to generate Airbus ALO problems of any size.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description