Method and System for solving problems with multiple conflicting objectives

Abstract

A computer implemented method for obtaining a set of Pareto-optimal solutions to a given multi-objective optimization problem with different objectives. The method generates a set of parameters on a classical computation component using the classical optimization algorithm. Based on the set of parameters, a quantum component generates a quantum state by executing a parameterized quantum circuit. A set of basis states is selected from the one quantum state The classical computation component calculates a set of objective values consisting of the objective values for each solution in the set of classical solutions, and a new estimate for the Pareto set of solutions and the corresponding Pareto front by the classical optimizer is generated. The internal state of the classical optimization algorithm is updated based on the current set of solutions and the Pareto set and their corresponding sets of objective values.

Description

CROSS-REFERENCE TO RELATED APPLICATION

This application claims the priority benefits of European application no. 23 181 719.8, filed on Jun. 27, 2023. The entirety of the above-mentioned patent application is hereby incorporated by reference herein and made a part of this specification.

The invention relates to a computer implemented method for obtaining a set of Pareto-optimal solutions to a given multi-objective optimization problem and a hybrid quantum-classical computer system implementing such method.

Quantum optimization is a very active field of research. The current generation of quantum computing hardware devices is still limited in sizes and subject to noise. The usefulness of these so-called noisy intermedia scale quantum (NISQ) devices for application problems is currently being tested and evaluated on many real-world application problems.

The current state of the art approaches for solving optimization problems with quantum computers are hybrid quantum-classical algorithms as described for example in “Noisy intermediate-scale quantum (NISQ) algorithms, Rev. Mod. Phys. 94, 015004”. Examples are variational quantum eigensolver (VQE) described in “The variational quantum eigensolver: a review of methods and best practices. Physics Reports, 986, 1-128”, quantum approximate optimization algorithm (QAOA) described in “A Quantum Approximate Optimization Algorithm, arXiv: 1411.4028 (2014)” or “From the quantum approximate optimization algorithm to a quantum alternating operator ansatz, Algorithms 12 (2019)” or the like as in “Quantum machine learning beyond kernel methods”, arXiv preprint arXiv: 2110.13162 (2021)”. All these approaches consider problems where the goal is to find the optimal solution(s) which minimize a predefined cost function. Technically, they target to solve single objective optimization tasks.

What is still missing is an equivalent quantum approach which targets the solutions of a multi-objective optimization problem.

There are a few approaches which actually utilize quantum hardware and consider multi-objective problems. In “Multiobjective variational quantum optimization for constrained problems: an application to Cash Management, arXiv: 2302.04196, https://doi.org/10.48550/arXiv.2302.04196”, a VQE quantum circuit and a coupled classical NSGAII algorithm is applied to it. From the quantum state prepared by the quantum apparatus, they calculate classically two objective functions: the cost function, which implements the actual goal when minimized and the probability to find feasible solutions in the final state. The target is not to map a full Pareto front, but arrive at one optimal solutions point, where the probability to find feasible solutions is very close to one and the cost function is minimized. In effect, this method is a constraint handling method, where a classical multi-objective optimization is used to efficiently explore the search space of feasible and infeasible solutions. No quantum aspect of the quantum computation is utilized to arrive at the solution of a multi-objective problem.

Despite the problems being attacked becoming more realistic in nature, one central feature of real-world problems is still not much considered: realistic problems often deal with multiple conflicting objectives such cost (which should be as low as possible), quality (which should be as high as possible), user satisfaction (which should be as large as possible), CO2 footprint (which should be as low as possible), . . . . These objectives are usually in a tradeoff relation where solutions which have very good values in one objective, c.g., high quality, are not very good in another objective, i.c., also have high cost. Therefore, the target of multi-objective optimization is to find a whole set of solutions (not only one) which give a good approximation to the so-called Parcto Front (PF), which is the best achievable set of trade-off set of solutions. Each solution in the PF is equal or better in at least one objective and at the same time worse in at least one other objective than any other solution in the PF.

“Multi-Objective Routing Optimization for 6G Communication Networks Using a Quantum Approximate Optimization Algorithm”, Sensors 22 (19), 7570 (2022) doi: 10.3390/s22197570 describes solving a multi-objective problem by iteratively solving many constrained single-objective problems. The quantum computer is only used to solve single-objective problems. Each single-objective problem optimizes only one objective while imposing constraints on the other objectives. It does not utilize the quantum computer to get a population of solutions directly.

Several papers exist which utilize so-called nondominated quantum optimization (NDQO) and modifications thereof like multi-objective decomposition quantum optimization (MODQO) or non-dominated quantum iterative optimization (NDQIO) described in “Quantum-Assisted Routing Optimization for Self-Organizing Networks”, IEEE Access 2, 614-632 (2014) doi: 10.1109/ACCESS.2014.2327596, “Non-Dominated Quantum Iterative Routing Optimization for Wireless Multihop Networks”, IEEE Access 3, 1704-1728 (2015) doi: 10.1109/ACCESS.2015.2478793 or “Quantum-Assisted Joint Multi-Objective Routing and Load Balancing for Socially-Aware Networks”, IEEE Access 4, 9993-10028 (2016) doi: 10.1109/ACCESS.2016.2629671. All these approaches target to find the Pareto optimal solution set to a multi-objective problem. They do so by utilizing a quantum computing apparatus. They specifically design a quantum circuit that evaluates Pareto dominance relations between individual solutions with the help of auxiliary qubits. Then they use variants/extensions of Grover search to search for novel good candidate solutions which they then can possibly insert into the current approximation of the Pareto front. These algorithms are specifically tailored to solve one single problem and use complex circuits to evaluate Pareto-dominance relations between sets of solutions. They cannot be applied to general multi-objective optimization problems directly, and in particular they cannot be executed on NISQ devices but require fully fault tolerant and error corrected quantum hardware not available in the near-to mid-term future.

However, there is still a need to directly address optimization of multiple conflicting objectives to enable an efficient calculation of a plurality of optimized solutions of the problem. The current state of art does not allow for handling multi-objective optimization problems natively on quantum hardware devices.

It is to be noted that in the following explanation, “classical computing component” may be a classical computer, and the “quantum computing component” or “quantum component” may be a universal or special-purpose purpose quantum computing device.

The current invention addresses the above mentioned problem to find solutions of realistic real world problems, specifically known from the energy management domain, on current and future quantum hardware devices employing a quantum-native multi-objective quantum optimization approach according to the method and computer system defined in the independent claims.

The dependent claims define advantageous aspects and features.

The invention encompasses the utilization of the quantum superposition, which is inherent to quantum states, as a resource to encode a complete population of solutions of a multi-objective optimization algorithm. The quantum circuit is designed in such a way as to increase the diversity in the population and also improve all objectives such that they give a good approximation to the Pareto front.

The method and apparatus for obtaining a preferably complete set of N_PF∈⁺ Parcto-optimal solutions, ={x_i}_i=1^NPF, to a given multi-objective optimization problem with M different objectives {right arrow over (C)}(x)={C_α}_α=1^M=(C₁(x), C₂(x), . . . , C_M(x))^T where x represents a classical solution of the multi-objective optimization problem, which is comprised of the following components and method steps:

• • a hybrid quantum-classical computer apparatus, with a quantum computing component and a classical component. The quantum computing component and the classical computing component may be realized as two separate entities, which are capable of exchanging information. Alternatively, they may be integrated into one device. The quantum computing component realizes a parametrized quantum circuit (PQC) by implementing a unitary quantum operation U(γ) which is parametrized by freely choosable parameters γ and which generates a N-qubit or qudit quantum state |ψ_γ by applying the unitary operator to a previously defined N-qubit or qudit initial state |ψ₀, i.e. |ψ_γ=U(γ)|ψ₀, which can be decomposed into the measurement basis of the quantum device, |ψ_γ=Σ_c=1^dHξ_γ,c|x_c with d_H the dimensionality of the Hilbert space, and each measurement basis state |x_c encodes one classical solution x_c,. The classical computing component is able to run a classical optimization algorithm by providing parameters values γ to the quantum computing component, initializing the execution of the PQC and the preparation of the state |ψ_γ, and then receiving a set of classical solutions from the quantum component as a result and calculating the cost functions for each solution. In essence, both components together are used to solve a multi-objective optimization problem by iteration. At each iteration t the following steps are performed:
  • a. generating a set of parameters γ_t on the classical computing component using the classical optimization algorithm and sending them to the quantum component,
  • b. on the quantum computing component, generating a quantum state |ψ_γt=U(γ_t)|ψ₀ by executing the PQC with the set of parameters γ_t, where the one quantum state |ψ_γt represents a (whole) set of different N_S solutions to the original problem, since it is naturally represented as a quantum mechanical superposition of many measurements basis states, i.e. |ψ_γt=Σ_c=1^dHξ_γt,c|x_c and each measurement basis state |x corresponds to one classical solution x of the original optimization problem,
  • c. using a basis-state selection method to select a set of N_S≥N_PF basis states {|x_i}_i=1^NS from the one quantum state |ψ_γt and sending the information of which basis states where selected, i.e., the set of classical solutions _t={x_i}_i=1^NS, to the classical computing component,
  • d. on the classical computing component calculating a set of objective values consisting of the objective values for each solution in the set of classical solutions, i.e. _t={{right arrow over (C)}(x_i)}_i=1^NS≡{right arrow over (C)}(_t),
  • e. generating a new estimate for the Pareto set of solutions, _t={x_i}_i=1^NPF, and the corresponding Pareto front, _t={{right arrow over (C)}(x_i)}_i=1^NPF, by the classical optimizer incorporating the set of solutions _t and the corresponding set of objectives values _t. According to a preferred embodiment, all those values from previous iterations, {_t′}_t′=1^t-1 and {_t′}_t′=1^t-1 which are retrieved from an archive are additionally used.
  • f. Based on the current set of solutions and the Pareto set and their corresponding sets of objective values, updating the internal state of the classical optimization algorithm and, in case that values from previous iterations shall be used in the generation of a new estimate, updating the archive =+{t, γ^t, _t, _t}
  • g. checking a convergence criterion as a stopping criterion, and if the stopping criterion is not satisfied, repeating the process from step a. onwards, until the stopping criterion is fulfilled.

The final results after the iterations have been completed is a good set of Pareto optimal solutions _Final={x_i}_i=1^NPF where the corresponding set of objective value vectors _Final={{right arrow over (C)}(x_i)}_i=1^NPF is an efficient and good approximation to the true Pareto front.

The invention uses quantum superpositions to encode a whole population of solutions into one quantum state. The quantum computing component is designed such that the quantum circuit layout achieves multi-objective evolution of Pareto-optimal solutions as a whole with a true quantum encoding.

According to a preferred embodiment, the quantum component realizes a parametrized quantum circuit (PQC) by implementing a unitary quantum operation U(γ) which is composed of N_L layers with different choosable sets of parameters as

$U\left(\gamma \right)=\prod _{l=1}^{N_L}{U}_l\left({\gamma }_l\right)$

and where each layer l comprises multiple mixing and phase operators with separately choosable parameter sets γ_l={y_X;l,γ_P;l}, and the unitary in each layer has the form

$U_l\left({\gamma }_l\right)=\prod _{\alpha =1}^M{e}^{-i{X}_{\alpha }\left({\gamma }_{X;l,\alpha }\right)}{e}^{-i{\gamma }_{H;l,\alpha }{\mathscr{H}}_{\alpha }}$

where the phase unitaries in each layer ˜ have freely choosable parameters γ_H;l,α∈ and are generated by the M cost Hamiltonians _α and α=1, . . . , M, where each cost Hamiltonian _α represents one classical objective function C_α, in the sense that the measurement basis states |x are simultaneous eigenstates of all problem Hamiltonians and the eigenvalues of the problem Hamiltonians are given by the value of the corresponding classical objective function of the corresponding classical solutions, x, i.e.

${\mathscr{H}}_{\alpha }❘x\rangle =C_{\alpha }\left(x\right)❘x\rangle$

and where the mixing unitaries in each layer, ˜e^{−iXα(γX;l,α)}, have freely choosable parameters γ_X;l,α and are generated by M parametrized mixing operators X_α(γ_X;l,α) with α=1, . . . , M, where the mixing operators do not commute with the cost Hamiltonians, i.e., [X_β(γ_X;l,β), H_α]≠0, ∀α, β, and therefore the mixing operators generate transitions between the measurement basis states, X_β(γ_X;l,β)βx_i=Σ_jζ_β,ij|x_j with ζ_β,ij∈.

According to another preferred embodiment, the quantum component realizes parametrized quantum circuit (PQC) by implementing a unitary quantum operation U(γ) which is composed of N_L layers with different choosable sets of parameters as

$U\left(\gamma \right)=\prod _{l=1}^{N_L}{U}_l\left({\gamma }_l\right)$

and where each layer l comprises multiple mixing and phase operators with separately choosable parameter sets γ_l∈{γ_X;l, γ_P;l},, and which has the form

$U_l\left({\gamma }_l\right)=\prod _{k=1}^{N_P}{e}^{{ }_{}-{\mathrm{iX}}_k\left({\gamma }_{X;l,k}\right)}{e}^{{ }_{}-i{\gamma }_{P;l,k}{P}_k}$

where the phase unitaries in each layer ˜e^{−iγP;l,kPk} have freely choosable parameters γ_P;l,k∈ and are generated by specific sums over the cost Hamiltonians _α

$P_k=\sum _{\alpha =1}^M{a}_{k,\alpha }{\mathscr{H}}_{\alpha }$

The coefficients in these sums a_K,α∈ may specify reference vectors {right arrow over (a)}_k which have to fulfill Σ_αa_k,α²={right arrow over (a)}_k^T·{right arrow over (a)}_k=1 and a_k,α≥0 for each k, and for each k, define a certain direction in the space of objective values and thus can either reflect one user preference and are determined in the beginning as part of the problem definition, or determined by the state of the art procedure for reference vectors as described in “A Reference Vector Guided Evolutionary Algorithm for Many-Objective Optimization”, IEEE Transactions on Evolutionary Computation 20(5), 773-791 (2016)doi 10.1109/TEVC.2016.2519378 or “.Adaptation of Reference Vectors for Evolutionary Many-objective Optimization of Problems with Irregular Pareto Fronts”, IEEE Congress on Evolutionary Computation (CEC), 1726-1733 (2019) doi 10.1109/CEC.2019.8790214, which can either be fixed reference vectors determined prior to iteration of the hybrid quantum-classical algorithms or adapted during the iterations in each iteration step t.

Each cost Hamiltonian _α represents one classical objective function C_α, in the sense that the measurement basis states |x are simultaneous eigenstates of all cost Hamiltonians and the eigenvalues of the problem Hamiltonians are given the value of the corresponding classical objective function of the corresponding classical solutions, x, i.e.

${\mathscr{H}}_{\alpha }❘x\rangle =C_{\alpha }\left(x\right)❘x\rangle$

and where the mixing unitaries in each layer, e^{−iXk(γX;l,k)}, have freely choosable parameters γ_X;l,k and are generated by M parametrized mixing operators X_α(γ_X;l,k) with k=1, . . . , N_P, where the mixing operators do not commute with the cost Hamiltonians, i.e., [X_k(γ_X;l,k),_α]≈0, ∀k, α, and therefore the mixing operators generate transitions between the measurement basis states, X_k(γ_X;l,k)|x_i=Σ_jζ_k,ij|x_j with ζ_k,ij∈. and where the N_S basis states which represent the set of solutions and which are encoded in the quantum state |ψ_γ as superpositions, are extracted at each iteration by estimating the full state tomography, i.e., by estimating the spectral weight ξ_γ,i with which each basis state |x_i contributes to the final state |ψ_γ, which is given by the absolute square of the overlap, i.e. ξ_γ,i=|x_i|ψ_γ|², and then selecting those N_S classical solutions states which have the largest spectral weights, and adding those to the set of selected solutions _t={x_i}_i=1^NS obtained in this iteration t.

According to another preferred embodiment the quantum component realizes parametrized quantum circuit (PQC) by implementing a unitary quantum operation U(γ) which is composed of N_L layers with different choosable sets of parameters as

$U\left(\gamma \right)=\prod _{l=1}^{N_L}{U}_l\left({\gamma }_l\right)$

and where each layer l comprises multiple mixing and phase operators with separately choosable parameter sets γ_l∈{γ_X;l,γ_V;l}, and which has the form

U _l(γ_l)=e^{−iXl(γX;l)}e^−iγV;lV

where the phase unitaries in each layer e^−iγV;lV have freely choosable parameters γ_V;l∈ and are generated by the operator V which estimates the hypervolume in objective space which is spanned by a population of states as it is encoded by a quantum superposition of states, i.e.

$V❘\psi \rangle =f\left({\left\{{\xi }_i\right\}}_{i=1}^{N_S}\right)\mathrm{HV}\left({\left\{\stackrel{\to }{C}\left(x_i\right)\right\}}_{i=1}^{N_s}\right)❘\psi \rangle$

where the function HV({{right arrow over (C)}(x_i)}_i^NS gives the hypervolume which is spanned by the objective value vectors {right arrow over (C)}(x_i)={C_α(x_i)}_α=1^M of the set of solutions {x_i}_i=1^NS and the function f({ξ_i}_i=1^NS) has function values between zero and one, i.e. 0≤f({ξ_i}_i=1^NS)≤1, and allows for a down-scaling of the hypervolume of a particular set of solutions in case the individual amplitudes of some of the contributing states are unfavorably small or the distribution of the amplitudes is too inhomogeneous, and where each cost Hamiltonian _α represents one classical objective function C_α, in the sense that the measurement basis states |x are simultaneous eigenstates of all cost Hamiltonians and the eigenvalues of the problem Hamiltonians are given the value of the corresponding classical objective function of the corresponding classical solutions, x, i.e.

_α |x =C _α(x)|x

and where the mixing unitaries in each layer e^{−iXl(γX;l)}, haves freely choosable parameters γ_X;l∈ and are generated by the parametrized mixing operator X_l(γ_X;l), where the mixing operators do not commute with the cost Hamiltonians, i.e., [X_l(γ_X;l), _α]≈0, ∀α, l, and therefore the mixing operators generate transitions between the measurement basis states, X_k(γ_X;l)|x_i=Σ_jζ_k,ij|x_j with ζ_k,ij∈

It is further preferred that the mixing operator in each layer is given by a sum over all x-angular momentum operators of each qudit variable and where the coefficient of each individual angular momentum operators is given by the freely choosable parameter γ_X;k∈, i.e.,

$X_k\left({\gamma }_{X;k}\right)=\sum _{i=1}^N{\gamma }_{X;k,i}{L}_{x;i},$

where L_x;i is the x-angular momentum operators of the i^th qubit or qudit with i=1, . . . , N, representing the i^th classical search variable and N the total number of qubits/qudits.

As an alternative embodiment, the mixing operator in each layer is given by a sum over all generalized Pauli X operators of each qudit variable as, for example, described in Y. Wang, Z. Hu, B. C. Sanders, and S. Kais, Qudits and High-Dimensional Quantum Computing, Frontiers in Physics 8 (2020), and where there is one coefficient for each physical qudit which are given by the freely choosable parameters γ_X;l,i∈, i.e.,

$X_k\left({\gamma }_{X;l,i}\right)=\sum _{i=1}^N{\gamma }_{X;l,i}{𝒳}_i$

where χ_i is the generalized Pauli X operator of the i^th qubit or qudit with i=1, . . . , N, and is given by χ_i=Σ_z=0^L-1|(z+1)mod Lz| and local qudit dimension L.

In a preferred embodiment the N_S basis states which represent the set of solutions and which are encoded in the quantum state |ψ_γ, as superpositions, are extracted at each iteration by sampling from the quantum state |ψ_γ, i.e. sampling x from the Born rule probability P(x, γ)=|x|ψ_γ|² which is realized by repeatedly measuring the quantum state, and adding for each measured state |x_m the corresponding classical solution x_m to the set of selected states, until N_S different states have been measured and added to the set of selected states, ={x_i}_i=1^NS.

The method and apparatus may further be applied to solve a multi-objective bi-directional electric vehicle charging problem where the target is to find a set ={x_i}_i=1^NPF of Pareto optimal configurations of discrete integer variables x=(₀, ₁, . . . , _Ntot), which encode a full charging and service schedule for a fleet of electric vehicles and where multiple objectives such as electricity cost for fulfilling all charging demands, customer satisfaction, maximizing utilization of renewable energies, minimizing electricity drawn from the grid, minimizing the peak electricity power consumption, etc are considered such that a good approximation of the Pareto front, is obtained.

The method and apparatus according to one of the above described configurations may be applied to obtain an efficient frontier of possible portfolio compositions as formulated by the Markowitz model of Modern portfolio theory by finding the compositions of a portfolio over time as specified by integer variables x={ω_n,s} which specify the investment amount into asses n at time s and where the efficient frontier is spanned by those portfolio compositions which have optimal tradeoffs between maximizing return and minimizing risk and thus pose a two-objective optimization problem.

The invention will now be further described with reference to the attached drawing, in which

FIG. 1 shows a schematic of a hybrid quantum-classical computer system illustrating the method steps according to a preferred implementation of the invention.

The proposed computer implemented method and hybrid quantum-classical computer system target to solve combinatorial integer multi-objective optimization problems with two or more objectives to be minimized. Typically, these objectives cannot be minimized simultaneously by some optimal solutions since they are in tradeoff relations to each other. Therefore, the target of the multi-objective optimization approach is to find a whole set of different solutions, called the Pareto set, where each solution realizes one optimal tradeoff between the objectives. The set of solutions therefore should map all possible tradeoff between the objectives, which is called the Pareto front.

The proposed inventive method and apparatus enables the treatment of such problems natively on current NISQ quantum hardware. The crucial aspect of the present invention is provided by representing a full population of candidate solutions to the optimization problem as one quantum state and extracting whole sets of Pareto-optimal solutions from this one quantum state.

According to the invention, the quantum circuits of the proposed computer implemented method and hybrid quantum-classical computer system increase the diversity of the population by utilization of all the cost function Hamiltonians corresponding to all classical objectives.

According to this description, the expression “classical” always refers to known computation as performed in the past using digital computer systems only.

The proposed computer implemented method and hybrid computer system comprises a classical computing component and a quantum computing component. Both components are configured to transfer and receive information from the respective other component such that a shared computation can be achieved in order to solve the multi objective integer optimization problem. A set of integer variables of the multi-objective integer optimization problem, x=(a₁, a₂, . . . , a_N), where a_i∈[0, . . . , L−1] are finite integers, are represented by the qubit (for L=2) or qudit (for L>2) degrees of freedom of the quantum hardware device, i.e.,

$❘x\rangle =❘a_1,a_n,...,a_N\rangle .$

Each of those states is a basis state in the measurement basis and fulfils the orthonormality conditions

$\langle x^{{ }_{}\prime }❘x\rangle =\langle a_1^{{ }_{}\prime },a_2^{{ }_{}\prime },\dots ,a_N^{{ }_{}\prime }❘a_1,a_2,\dots ,a_N\rangle ={\delta }_{a_1^{{ }_{}\prime },a_1}\dots {\delta }_{a_N^{{ }_{}\prime },a_N}$

where δ_b,a is the Kronecker delta,

${\delta }_{b,a}=\{\begin{array}{cc}1& \mathrm{if}a=b\\ 0& \mathrm{if}a\ne b\end{array}.$

This encoding means that all possible solutions to the original multi objective integer optimization problem are represented by the measurement basis states of the quantum computing component. One advantage of using a quantum computer for approaching this problem is given by the fact that a general quantum state is naturally a superposition of all basis states, i.e. a superposition of all possible solutions to the problem,

$❘\varphi \rangle =\sum _{c=1}^{d_H}{\xi }_c❘x_c\rangle$

with complex coefficient ξ_c∈ and where Σ_C|ξ_c|²=1 and where c=1, . . . , d_H enumerates all basis states of the Hilbert space of dimensions d_H=L^N for N variables (i.e. qubits or qudits) with L states per variable (i.e. qubits or qudits).

The invention utilizes a parametrized quantum circuit on the quantum computing component, which employs unitary gate operations that are applied to a predefined fixed initial state |ψ₀. Such application of the unitary gate operations leads to the final quantum state

|ψ=U(γ)|ψ₀.

According to the invention, this state encodes a full population of N_S solutions candidates of the original problem,

$❘\psi \rangle =\sum _{i=1}^{N_s}{\xi }_i❘x_i\rangle +\sum _{j=N_S+1}^{d_H}{\delta }_j❘x_j\rangle$

where the first sum are those basis states with large weights and which represents the full population, ={x_i}_i=1^NS, and the second sum is negligibly small, i.e. Σ_i=1^NS|ξ_i|²=1−ϵ and Σ_j=NS+1^dH|δ_j|²=ϵ<<1, respectively, such that the final state truly encodes the solution population,

$❘\psi \rangle \approx \sum _{i=1}^{N_s}{\xi }_i❘x_i\rangle$

A state selection mechanism, such as full state tomography, classical shadows or measurement-based sampling allows the full set of solutions to be extracted from the quantum state,

$❘\psi \rangle 𝕊={\left\{x_m\right\}}_{m=1}^{N_S}.$

Suitable techniques are described in “Efficient quantum state tomography”, Nat Commun 1, 149 (2010) doi 10.1038/ncomms114, “Predicting many properties of a quantum system from very few measurements”, Nat. Phys. 16, 1050-1057 (2020) doi 10.1038/s41567-020-0932-7 or “Quantum Computation and Quantum Information”, Cambridge: Cambridge University Press. (2010) doi 10.1017/CBO9780511976667. The unitaries of the current invention are specifically designed to increase the contributions to the superposition from different states, which provide the optimal trade-off solutions. That means, the population of candidate solutions represented by one quantum state will provide a good approximation to the Pareto set of solutions ,

$❘\psi \rangle =\sum _{i=1}^{N_{\mathrm{PF}}}{\xi }_i❘x_i\rangle \iff \mathrm{\mathbb{P}𝕊}={\left\{x_i\right\}}_{i=1}^{N_{\mathrm{PF}}}.$

In the proposed invention, the set of solutions will be sent from the quantum computing component to the classical computing component, where for each solution the vector of objective values will be calculated which generates the Pareto front,

$\mathrm{\mathbb{P}𝔽}={\left\{\stackrel{\to ︀}{C}\left(x_i\right)\right\}}_{i=1}^{N_{\mathrm{PF}}}.$

Here, the vector {right arrow over (C)}(x) represents the classical M objective functions of the original problem,

$\stackrel{\to ︀}{C}\left(x\right)={\left\{C_{\alpha }\right\}}_{\alpha =1}^M={\left(C_1\left(x\right),C_2\left(x\right),\dots ,C_M\left(x\right)\right)}^T.$

The final result after executing the proposed method using the inventive hybrid quantum-classical computing system is the Pareto set of solutions , which approximate the Pareto front ={{right arrow over (C)}(x_i)}_i=1^NPF. In that way, the user can inspect the Pareto optimal solutions and then select the best tradeoff solutions for their application in mind. The proposed method and apparatus utilizing the PQC is effective in producing a good approximation to the Pareto front since all possible solutions are typically present in the superposition of the initial state |ψ₀ of the quantum computing component. The proposed phase operators in the PQC do assign different phases to each individual basis state |ψ. In that way, the classical optimizer can tune the parameters of the PQC to give the largest diversity and include mostly non-dominated solutions in the quantum state, which maximize hypervolume indicator which implies a good approximation to the true Pareto front.

Hereinafter, the invention will be described for two different implementation examples, which are typical for multi objective problems:

The operation of a fleet of electric vehicles (EVs) constitutes an archetypical optimization problem with several challenges. The central aspect of this operation is to provide a service to the customer such that for all vehicles, the requests of the customers can be fulfilled while at the same time reducing the electricity cost for charging all EVs as much as possible.

In one formulation of the problem, the target is to find the optimal charging schedule of a fleet of N_EV electric vehicles (EVs) over N_T time slots. The search variable x defines the charging schedule of the complete EV fleet. A schedule consists of N_EV·N_T discrete integer variables,

x=(_1,1, _2,1, . . . , _EV,1, _2,1, . . . , _n,s, . . . , _NEV,NT)

where each integer variable can have only one of L values, _n,s∈[0,1, . . . , L−1]. The value of the integer _n,s determines the charging power level of vehicle n in time slot s, i.e. p (s, n) =p_0n,s where the p₀ is the conversion constant and thus defines the spacing between the charging levels. There are four common objectives, which should all be simultaneously minimized. These are expressed by the following four objective functions:

$\begin{array}{c}{C}_1\left(x\right)=\frac{1}{N_1}\sum _{s=1}^{N_T}𝒫\left(s\right)\left(\sum _{n=1}^{N_{\mathrm{EV}}}{p}_0{\ell }_{n,s}-\mathrm{PV}\left(s\right)\right)\left(\mathrm{OBJ}1,\mathrm{total}\mathrm{cost}\right)\\ C_2\left(x\right)=\frac{1}{N_2}\sum _{n=1}^{N_{\mathrm{EV}}}{\left(\sum _{s=1}^{N_T}{p}_0{\ell }_{n,s}-E_{\mathrm{tot},n}\right)}^2\left(\mathrm{OBJ}2,\mathrm{fairness}\right),\\ C_3\left(x\right)=\frac{1}{N_3}\sum _{s=1}^{N_T}{\left(\sum _{n=1}^{N_{\mathrm{EV}}}{p}_0{\ell }_{n,s}-\mathrm{PV}\left(s\right)\right)}^2\left(\mathrm{OBJ}3,\mathrm{PV}\mathrm{usage}\right),\\ C_4\left(x\right)=\frac{1}{N_4}\sum _{s=1}^{N_T}{\left(\sum _{n=1}^{N_{\mathrm{EV}}}{p}_0{\ell }_{n,s}\right)}^2\left(\mathrm{OBJ}4,\mathrm{power}\mathrm{variations}\right).\end{array}$

with the normalization constants N₁=N_TN_EVmax/p₀, N₂=N_T²N_EV(L−1)/p₀², N₃=N_TN_EV²/p₀², and N₄=N_TN_EV²/p₀² to ensure that the numerical values of the objectives are of approximately the same order of magnitude.

The first objective (OBJ1) represents to total electricity cost for the complete charging schedule. PV(s) denotes the surplus energy produced by a photo-voltaic (PV) system at time slot s which is available for charging the EVs, and the function (s) gives the price for buying or selling energy from or to the grid at time slot s. _max denotes the maximal unit-price for electricity. When the total electricity demand is larger than the PV surplus, energy needs to be bought from the grid, while when is smaller the term in brackets reverses sign end energy is sold to the grid leading to a negative cost, i.e. income.

The second objective (OBJ2) represents a measure for fairness between the different

EVs/customers by minimizing the difference between the desired energy of each vehicle, E_tot,n, and the energy provided to each EV/customer. This is achieved since this term is proportional to the variance of the difference between the charging energy provided to each EV n and its desired value, E_tot,n, which strongly favors equal distribution.

The third objective (OBJ3) tries to only use the PV surplus energy for charging the EVs at each time slot s, and consequently the objective function tries to minimize the difference between used surplus PV energy and current total charging energy at each time slot. Minimizing this objective is a means to make the service operation independent of the electricity grid.

The fourth objective (OBJ4) targets the equal distribution of the charging power over the time slots, and in particular tries to reduce peaks in the charging power at all time slots. The reasoning behind this is two-fold: (i) for company settings, the peak electricity power over a billing period has a direct influence on the charged electricity bill and therefore strong demand peaks need to be avoided (ii) strong variations in the electricity demand generally put a lot of pressure on the energy grid provider which needs to balance these out. Both aspects are met by reducing the overall variation over time of the total energy demand, i.e. the variance, as formulated in this objective.

The problem characterized by the objectives (OBJ1)-(OBJ4) of the above described example defines a combinatorial integer multi-objective optimization problem with M=4 objectives, which all need to be minimized. However, the present invention is not limited to four objectives and less or more objectives may define the multi objective optimization problem. These objectives typically cannot be minimized simultaneously by some optimal solutions since they are in tradeoff relations to each other. Therefore, the target of the multi-objective optimization approach according to the invention is to provide a computationally efficient way to find a whole set of different solutions, called the Pareto set, where each solution realizes one optimal tradeoff between the objectives. The set of solutions therefore should map all possible tradeoff between the objectives which is called the Pareto front.

For the EV charging problem the discrete search variables are given by the variables x=(_1,1,_2,1, . . . , _n,1,_2,1, . . . ,_n,s, . . . ,_NEV,NT) with _n,s∈[0,1, . . . , L−1] which encode the full charging schedule of all EVs. The classical cost functions, i.e. the objectives to be minimized are expressed with those variables, as is exemplarily done in (OBJ1)-(OBJ4).

In the EV charging problem the encoding into the quantum formulation is as follows: each qudit of the quantum computing component represents one charging power level _n,s∈[0, . . . , L−1] of one Ev n∈[1, . . . , N_EV] in one time slot s∈[1, . . . , N_T]. The quantum states of the measurement basis are given by

$\begin{array}{cc}❘x\rangle =❘{\ell }_{1,1},{\ell }_{2,1},\dots ,{\ell }_{N_{\mathrm{EV}},1},{\ell }_{2,1},\dots ,{\ell }_{n,s},\dots ,{\ell }_{N_{\mathrm{EV}},N_T}\rangle & \left(\mathrm{ENC}\right)\end{array}$

with _n,s∈[0,1, . . . , L−1] which represent all the different possibilities for the charging schedule of the whole EV fleet. This encoding implies that all possible solutions to the EV charging problem are represented by the measurement basis states of the quantum circuit. One advantage of using a quantum computer for approaching this problem is given by the fact that a general quantum state in this basis is naturally a superposition of all basis states, i.e. a superposition of all possible solutions to the problem,

$❘\varphi \rangle =\sum _x{c}_x❘x\rangle$

with complex coefficient c_x∈ and where Σ_x|c_x|²=1.

The classical M objective functions C₁(x) to C₄(x) as defined above are represented by M cost Hamiltonians. Their functional form is given by a quadratic polynomial in the search variables. This directly allows to formulate the corresponding quantum Hamiltonians for each objective by taking the classical objective functions C₁(x) to C₄(x) and replacing the search variables _n,s which specify the charging level of EV n in times slot s with the corresponding z—angular momentum operator L_z;n,s of the qudit encoding the corresponding variable, i.e. by using the replacement rule

${\ell }_{n,s}\to L_{z;n,s}+\frac{L-1}{2}.$

This leads to cost Hamiltonians ₁ to ₄ which are also quadratic polynomials in the angular momentum operators, which therefore can be directly implemented with the elementary gate- or pulse-operators on current quantum hardware without further pre-processing or approximation. The quadratic form is therefore the starting point for current state of the art quantum approaches, as it can be found for example in “Encoding-Independent Optimization Problem Formulation for Quantum Computing”, arXiv: 2302.03711 (2023) or “Noisy intermediate-scale quantum (NISQ) algorithms”, Rev. Mod. Phys. 94, 015004.

The measurement basis states are eigenstates to the z—angular momentum operators L_z;n,s and the eigenvalues are the charging levels, i.e.

$\left(L_{z;n,s}+\frac{L-1}{2}\right)❘{\ell }_{1,1},\dots ,{\ell }_{n,s},\dots ,{\ell }_{N_{\mathrm{EV}},N_T}\rangle ={\ell }_{n,s}❘{\ell }_{1,1},\dots ,{\ell }_{n,s},\dots ,{\ell }_{N_{\mathrm{EV}},N_T}\rangle .$

With this, the explicit form of the corresponding cost Hamiltonians representing the objective functions are

$\begin{array}{c}{\mathscr{H}}_1=\frac{1}{N_1}\sum _{s=1}^{N_T}𝒫\left(s\right)\left(\sum _{n=1}^{N_{\mathrm{EV}}}{p}_0·\left(L_{z;n,s}+\frac{L-1}{2}\right)-\mathrm{PV}\left(s\right)\right)\left(\mathrm{HAM}1,\mathrm{total}\mathrm{cost}\right)\\ {\mathscr{H}}_2=\frac{1}{N_2}\sum _{n=1}^{N_{\mathrm{EV}}}{\left(\sum _{s=1}^{N_T}{p}_0·\left(L_{z;n,s}+\frac{L-1}{2}\right)-E_{\mathrm{tot},n}\right)}^2\left(\mathrm{HAM}2,\mathrm{fairness}\right),\\ {\mathscr{H}}_3=\frac{1}{N_3}\sum _{s=1}^{N_T}{\left(\mathrm{PV}\left(s\right)-\sum _{n=1}^{N_{\mathrm{EV}}}{p}_0·\left(L_{z;n,s}+\frac{L-1}{2}\right)\right)}^2\left(\mathrm{HAM}3,\mathrm{PV}\mathrm{usage}\right),\\ {\mathscr{H}}_4=\frac{1}{N_4}\sum _{s=1}^{N_T}{\left(\sum _{n=1}^{N_{\mathrm{EV}}}{p}_0·\left(L_{z;n,s}+\frac{L-1}{2}\right)\right)}^2\left(\mathrm{HAM}4,\mathrm{power}\mathrm{variations}\right).\end{array}$

The measurement basis states are also simultaneous eigenstates of these cost Hamiltonians (since the latter are polynomials in the L_z;n,s), and the eigenvalues are given by the value of the objective functions of the corresponding configuration,

${\mathscr{H}}_{\alpha }❘{\ell }_{1,1},...,{\ell }_{n,s},...,{\ell }_{N_{\mathrm{EV}},N_T}\rangle =C_{\alpha }\left({\ell }_{1,1},...,{\ell }_{n,s},...,{\ell }_{N_{\mathrm{EV}},N_T}\right)❘{\ell }_{1,1},...,{\ell }_{n,s},...,{\ell }_{N_{\mathrm{EV}},N_T}\rangle .$

The quantum computing component, which executes the parametrized quantum circuit (PQC), prepares the state,

$\begin{array}{cc}❘\psi \rangle =\prod _{l=1}^{N_L}{U}_l\left({\gamma }_l\right)❘{\psi }_0\rangle & \left(PQC\right)\end{array}$

which consists of N_L different unitary operations U_l called layers and where each unitary operator is freely parametrizable with a set of parameters γ_l (l=1, . . . , N_L). The initial state of the circuit is given by the equal superposition of all possible states,

$❘{\psi }_0\rangle =\frac{1}{\sqrt{d_H}}\sum _x❘x\rangle ,$

where d_H=L^NEVNT is the dimensionality of the Hilbert space for N_EV EVs, N_T time slots, and where each time slot can be in one of L charging levels, encoded by N_EV·N_T qudits with local dimensionality L. In the current example, the unitary operation in each layer l has the form

$U_l\left({\gamma }_l\right)=\prod _{k=1}^{N_P}{e}^{-iX\left({\gamma }_{X;l,k}\right)}{e}^{-i{\gamma }_{P;l,k}{P}_k}$

which is a product over N_p different terms and each product term consists of a so-called mixing term, ˜e^{−iX(γx;l,k)} and a phase operator ˜e^{−iγP;l,kPk}.

The mixing unitary has N_PN_EVN_T freely choosable parameters γ_X;l,k,n,s∈ and it is generated by the Hermitian mixing operator X. The mixing operator in this example is given by the sum over all x-angular momentum operators of each physical qudit, i.e.,

$\begin{array}{cc}X\left({\gamma }_{X;l,k}\right)=\sum _{s=1}^{N_T}\sum _{n=1}^{N_{\mathrm{EV}}}{\gamma }_{X;l,k,n,s}{L}_{x;n,s},& \left(\mathrm{MIXEX}\right)\end{array}$

where L_x;n,s is the usual x-angular momentum operator of the physical qudit encoding the charging power at time slot s of EV n. The action of L_x;n,s is given by the standard action of angular momentum operators,

$L_{x;n,s}❘{\ell }_{1,1},\dots ,{\ell }_{n,s},\dots ,{\ell }_{N_{\mathrm{EV}},N_T}\rangle =\frac{1}{2}\sqrt{\left(L-{\ell }_{n,s}-2\right)\left({\ell }_{n,s}+2\right)}❘{\ell }_{1,1},\dots ,{\ell }_{n,s}+1,\dots ,{\ell }_{N_{\mathrm{EV}},N_T}\rangle +\frac{1}{2}\sqrt{\left(L-{\ell }_{n,s}\right){\ell }_{n,s}}|{\ell }_{1,1},\dots ,{\ell }_{n,s}-1,\dots ,{\ell }_{N_{\mathrm{EV}},N_T}\rangle$

and as such generates transitions between and creates superpositions of basis states, i.e. different solutions.

Each of the N_P phase terms in each layer comes with one freely chooseable parameter, γ_P;l,k∈ and each term is generated by a Hermitian operator which is one specific sum over the four cost Hamiltonians of (HAM1)-(HAM4), i.e.

$\begin{array}{cc}{P}_k=\sum _{\alpha =1}^4{a}_{k,\alpha }{\mathscr{H}}_{\alpha ,},& \left(\mathrm{PHASEEX}\right)\end{array}$

The coefficients a_k,α≥0 define reference vectors {right arrow over (a)}_k which need to fulfill the condition Σ_α=1⁴a_k,α²={right arrow over (a)}_k^T·{right arrow over (a)}_k=1, i.e., the vectors {right arrow over (a)}_k be unit vectors and thus each {right arrow over (a)}_k defines a direction in the objective space. These are the quantum analog of the reference vectors in classical multi-objective optimization as described in “A Reference Vector Guided Evolutionary Algorithm for Many-Objective Optimization”, IEEE Transactions on Evolutionary Computation 20(5), 773-791 (2016) doi 10.1109/TEVC.2016.2519378. In the current example it is assumed that N_P=5 (i.e. k=1, . . . ,5) reference vectors are used and the reference vectors {right arrow over (a)}_k are determined by one of the standard approaches from literature as described in, for example, “A Reference Vector Guided Evolutionary Algorithm for Many-Objective Optimization”, IEEE Transactions on Evolutionary Computation 20(5), 773-791 (2016) doi 10.1109/TEVC.2016.2519378, or “.Adaptation of Reference Vectors for Evolutionary Many-objective Optimization of Problems with Irregular Pareto Fronts”, IEEE Congress on Evolutionary Computation (CEC), 1726-1733 (2019) doi 10.1109/CEC.2019.8790214.

In a simple but representative case the quantum circuit has only one layer, N_L=1, in which case the quantum circuit prepares the following state,

$\begin{array}{cc}❘\psi \left({\left\{{\gamma }_{X;k},{\gamma }_{P;k}\right\}}_{k=1}^{N_P=5}\right)\rangle =\prod _{k=1}^5{e}^{-i{\sum }_{s=1}^{N_T}{\sum }_{n=1}^{N_{\mathrm{EV}}}{\gamma }_{X;k,n,s}{L}_{x;n,s}}{e}^{-i{\gamma }_{P;k}{\sum }_{\alpha =1}^4{a}_{k,\alpha }{\mathscr{H}}_{\alpha }}❘{\psi }_0\rangle & \left(\mathrm{PQCEX}\right)\end{array}$

with the reference vectors a_k,α and the cost Hamiltonians of (HAM1)-(HAM4). This parametrized quantum state depends on 5(N_TN_EV+1) parameters {{γ_X;k,n,s}_n,s^NEV,NT, γ_P;k}_k=1^NP=5 which need to be optimally determined by a classical optimization algorithm.

The claimed method executed on the claimed hybrid quantum-classical computer system comprises the following steps. A classical optimization algorithm such as a stochastic gradient descent (SGD) is initiated on the classical computing component. The search parameters of this classical optimizer are given by the 10 parameters of the PQC, i.e., γ={{γ_X;k,n,s}_n,s^NEV,NT, γ_P;k}_k=1^NP=5∈^D with D=5(N_TN_EV+1).

Additionally, an archive may be initiated on the classical computing component where search parameters, candidate solutions and corresponding values of the objective functions will be stored. Further, the quantum apparatus is configured such that the measurement basis represents classical solutions as depicted in (ENC) and the quantum circuit produces the parametrized quantum state (PQCEX) with the corresponding gate operations (MIXEX) and (PQCEX) and (HAM1)-(HAM4). Then, an iteration counter is initialized with t=0 and the following steps are iterated:

• • i) Based on its internal state, the classical optimizer algorithm generates a set of search parameter values,

$\mathrm{𝒞𝒪𝒜}\to {\gamma }^t={\left\{{\left\{{\gamma }_{X;k.n.s}^t\right\}}_{n,s}^{N_{\mathrm{EV}},N_T},{\gamma }_{P;k}^t\right\}}_{k=1}^{N_P=5}$

and the classical computing components sends these parameter values to the quantum apparatus.

• • ii) The quantum apparatus receives the parameter set γ^t={{γ_X;k,n,s^t}_n,s^NEV,NT, γ_P;k^t}_k=1^NP=5 from the classical computing component and adjusts its unitary gate operators such that it produces the quantum state specified by (PQCEX), i.e.

$❘\psi \left({\gamma }^t\right)\rangle =\prod _{k=1}^{N_P=5}{e}^{-i{\sum }_{s=1}^{N_T}{\sum }_{n=1}^{N_{\mathrm{EV}}}{\gamma }_{X;k.n.s}^t{L}_{x;n,s}}{e}^{-i{\gamma }_{P;k}^t{\sum }_{\alpha =1}^4{a}_{k,\alpha }{\mathscr{H}}_{\alpha }}❘{\psi }_0\rangle ,$

when it is executed.

Then the quantum apparatus executes the quantum circuit. The one resulting quantum state |ψ(γ^t) then encodes a whole population of N_S candidate solutions as superposition,

$❘\psi \left({\gamma }^t\right)\rangle =\sum _{i=1}^{N_s}{\xi }_i\left({\gamma }^t\right)❘x_i\rangle +\sum _{j=N_S+1}^{d_H}{\delta }_j\left({\gamma }^t\right)❘x_j\rangle$

where the first sum are those basis states with large weights and which represents the full population, ={x_i}_i=1^NS, and the second sum is negligibly small, i.e. Σ_i=1^NS|ξ_i(γ^t)|²=1−ϵ and Σ_j=NS+1^dH|δ_j(γ^t)|²=ϵ with ϵ much smaller than 1, respectively, such that the final state truly encodes the solution population,

$❘\psi \left({\gamma }^t\right)\rangle \approx \sum _{i=1}^{N_s}{\xi }_i\left({\gamma }^t\right)❘x_i\rangle .$

The set of N_S solutions _t at iteration t, is extracted from one quantum state as those Ns basis states which have the largest amplitudes.

One way to extract those solutions from the quantum state practically can be through repeated sampling. There, the quantum state is prepared and measured with the same parameters γ^t a plurality of times (i.e., N_measure>N_S), where each measurement results in one basis state, i.e., one solution,

$|\psi \left({\gamma }^t\right)\rangle \stackrel{\mathrm{measurement}\#m}{⟹}|x_m\rangle \stackrel{\mathrm{readout}}{\Rightarrow }{x}_m$

and collecting the measurement results in a list, from which the N_S most frequent measurement results are then selected. For large enough number of measurements, this sampling procedure indeed selects those basis states with the largest amplitudes, as each individual measurement selects a basis state according to the standard Born-rule probability,

Probability of measuring |x: P(x, γ^t)=|x|ψ(γ^t)|².

Another way could be to perform a full state tomography protocol or classical shadow tomography (see e.g. “Efficient quantum state tomography” Nat Commun 1, 149 (2010), doi 10.1038/ncomms1147, or “Predicting many properties of a quantum system from very few measurements”, Nat. Phys. 16, 1050-1057 (2020). https://doi.org/10.1038/s41567-020-0932-7) and explicitly selecting the basis states with the largest probabilities.

In any case, the procedure is such that the with the present invention, the quantum computing component extracts a whole population _t={x_m}_m=1^NS, of solutions from the one quantum state |ψ(γ^t) and then sends it to the classical computing apparatus for further processing.

• • iii) The classical computing component receives the population of solutions _t={x_m}_m=1^NS from the quantum computing component and stores it for further processing in an archive . Then, all objective function values from equations (OBJ1) to (OBJ4) are computed for each of the solutions in the population, which results in a set of objective value vectors _t={{right arrow over (C)}(x_m)∀x_m∈_t}. From this set of objective value vectors, the set of non-dominated solutions are selected, _t={{right arrow over (C)}(x_i)}_i=1^NPF. This represents the approximation of the Pareto front at iteration t and the corresponding set of solutions is an estimate of the Pareto set _t={x_i}_i=1^NPF.

Form this set of nondominated solutions, the hypervolume indicator HV(_t) is calculated, which in turn is used to calculate the fitness value f associated with the parameters of the PQC,

$f_t\equiv f\left({\gamma }^t\right)=-\mathrm{HV}\left({\mathrm{\mathbb{P}𝔽}}_t\right).$

The fitness value f_t is then fed back into the classical optimization algorithm_

←f(γ^t)

and the optimizer updates its internal state accordingly.

Additionally, the search parameters γ^t, the corresponding fitness value f(γ^t), the current estimate of the Pareto front _t, and the corresponding Pareto set of candidate solutions _t, are added to the archive,

→+{t, γ^t, f(γ^t), _t, _t}

• • iv) Thereafter, a convergence criterion is checked on the classical computing component. This can be, for example, whether the maximally allowed number of cost function evaluations has been reached, or if the improvement of the fitness value f(γ^t) over the last couple of iterations was not enough compared to a given threshold.
    • a. If the convergence criterion is met the optimization is stopped and as the final result the final Pareto front and corresponding Pareto set are obtained by selecting all the non-dominated solutions from all approximations of the Pareto front of each iteration _Final=NonDominatedSolutions({_t′}_t′=1^t).
    • b. If the convergence criterion is not met, the process is repeated from step i). above.

The final result is the Pareto set of solutions _Fina={x_i}_i=1^NPF which approximate the Pareto from _Final={{right arrow over (C)}(x_i)}_i=1^NPF, which is output for further evaluation or control of a charging system.

The proposed method and hybrid computer system utilizing the PQC is effective in producing a good approximation to the Pareto front since all possible solutions are typically present in the superposition of the initial state of the quantum apparatus. The phase operators in the PQC do assign different phases to each individual basis state based on the values of the reference vectors used in the phase operators of Eq. (PQCEX). In that way, the classical optimizer can tune the parameters of the PQC to give the largest diversity and include mostly non-dominated solutions in the quantum state which maximize hypervolume indicator.

The set of Pareto optimal solutions serves as a basis for implementing an efficient control of the EV charging and scheduling service. A solution selection routine can be implemented, where the best fitting solutions are automatically determined from the whole set of Pareto optimal solutions. The natural solution selection routine is to utilize state of the art knee point detection methods as, for example, described in “Performance Analysis on Knee Point Selection Methods for Multi-Objective Sparse Optimization Problems”, IEEE Congress on Evolutionary Computation, 1-8 (2018) doi 10.1109/CEC.2018.8477915 or “A Multi-objective Evolutionary Algorithm for “Finding Knee Regions Using Two Localized Dominance Relationships”, IEEE Transactions on Evolutionary Computation 25(1), 145-158 (2021) doi 10.1109/TEVC.2020.3008877. Knee-points naturally provide the preferred best tradeoff solutions.

The second application example regards a Markowitz portfolio optimization. The target of Markowitz portfolio optimization is to manage a portfolio of financial assets over a period of time such that the expected return is maximized, while at the same time the risk is minimized (see e.g., “Computational aspects of alternative portfolio selection models in the presence of discrete asset choice constraints”, Quantitative finance 1(5), 489-501 (2001), “Solving the Optimal Trading Trajectory Problem Using a Quantum Annealer”, IEEE Journal of Selected Topics in Signal Processing 10(6), 1053-1060, (2016) doi 10.1109/JSTSP.2016.2574703, or “Dynamic portfolio optimization with real datasets using quantum processors and quantum-inspired tensor networks” Phys. Rev. Research 4, 013006 (2022) for recent introductions). As these two objectives are usually in a tradeoff relation, it is desired to map the efficient frontier, i.e., the Pareto front of best possible tradeoff solutions.

Technically the problem consists of determining the compositions of a portfolio of financial assets as specified by the share ω_n,s of asset n, with n=1, . . . , N_Asset at time s with s=1, . . . , N_T for a fixed time window. Usually rebalancing, i.e., changing the portfolio composition is done a fixed predefined times which leads to these discretized time steps. Additionally, it is common for funds to trade assets in discrete junks, which leads to the shares ω_n,s to be representable by integer variables, i.e. ω_n,s∈[0, . . . , L−1].

The two objectives of this problem are formulated as

$\begin{array}{cc}{C}_1\left(\omega \right)=\sum _{s=1}^{N_T}\sum _{n=1}^{N_{\mathrm{Asset}}}\left[{\mu }_{n,s}{\omega }_{n,s}-v_{n,s}❘{\omega }_{n,s+1}-{\omega }_{n,s}❘\right]& \left(\mathrm{OBJ}1,\mathrm{Expected}\mathrm{return}\right)\end{array}$ $\begin{array}{cc}{C}_2\left(\omega \right)=\sum _{s=1}^{N_T}\sum _{n=1}^{N_{\mathrm{Asset}}}{\omega }_{n,s}{\sum }_{\mathrm{nm},s}{\omega }_{m,s}& \left(\mathrm{OBJ}2,\mathrm{Risk}\right)\end{array}$

Here μ_n,s is the forecasted return on asset n at time s and v_n,s is the transaction cost when changing the share of asset n from time step s−1 to time step s in the portfolio.

The first term thus gives the accumulated overall expected return for the complete portfolio over the considered time horizon.

In the second objective Σ_nm,s denotes the covariance of assets n and m at time s and the objective represents the covariance of the expected portfolio return summed over all times. This gives a measure for risk since (a) large covariance implies that the true return can differ largely form the expected returns, while small covariance implies the true return is similar to the expected return, and (b) large covariance between two assets implies correlations between the assets which in-turn reduces the robustness of the portfolio and thus makes the portfolio more risky.

The two-objective portfolio optimization problem is approached with the proposed method and hybrid quantum-classical computing system by representing the integer variables describing the shares in the portfolios, ω_n,s, as the measurement basis states of the quantum apparatus,

$\begin{array}{cc}❘x\rangle =❘{\omega }_{1,1},{\omega }_{2,1},...{\omega }_{N_{\mathrm{Assets}},1},{\omega }_{1,2},...,{\omega }_{n,s},...,{\omega }_{N_{\mathrm{Assets}},N_T}\rangle & \left(\mathrm{ENC}\right)\end{array}$

and formulating the two cost Hamiltonians in terms of the z-angular momentum operator L_z;n,s of the qudit encoding the corresponding variable, i.e. be using the replacement rule

${\omega }_{n,s}\to L_{z;n,s}+\frac{L-1}{2}.$

The measurement basis states are eigenstates to these operators and the eigenvalues are the shares of the asset in the portfolio, i.e.

$\left(L_{z;n,s}+\frac{L-1}{2}\right)❘{\omega }_{1,1},{\omega }_{2,1},...{\omega }_{N_{\mathrm{Assets}},1},{\omega }_{1,2},...,{\omega }_{n,s},...,{\omega }_{N_{\mathrm{Assets}},N_T}\rangle ={\omega }_{n,s}❘{\omega }_{1,1},{\omega }_{2,1},...{\omega }_{N_{\mathrm{Assets}},1},{\omega }_{1,2},...,{\omega }_{n,s},...,{\omega }_{N_{\mathrm{Assets}},N_T}\rangle .$

With this, the cost Hamiltonians are written as

$\begin{array}{cc}{\mathscr{H}}_1=\sum _{s=1}^{N_T}\sum _{n=1}^{N_{\mathrm{Asset}}}\left[{\mu }_{n,s}\left(L_{z;n,s}+\frac{L-1}{2}\right)-v_{n,s}❘L_{z;n,s+1}-L_{z;n,s}❘\right]& \left(\mathrm{HAM}1\right)\end{array}$ $\begin{array}{cc}{\mathscr{H}}_2=\sum _{s=1}^{N_T}\sum _{n=1}^{N_{\mathrm{Asset}}}\left(L_{z;n,s}+\frac{L-1}{2}\right){\sum }_{\mathrm{nm},s}\left(L_{z;m,s}+\frac{L-1}{2}\right).& \left(\mathrm{HAM}2\right)\end{array}$

The measurement basis states are also simultaneous eigenstates of these cost Hamiltonians and the eigenvalues are given by the value of the objective functions of the corresponding configuration,

${\mathscr{H}}_{\alpha }❘{\omega }_{1,1},...,{\omega }_{n,s},...,{\omega }_{N_{\mathrm{Assets}},N_T}\rangle =C_{\alpha }\left({\omega }_{1,1},...,{\omega }_{n,s},...,{\omega }_{N_{\mathrm{Assets}},N_T}\right)❘{\omega }_{1,1},...,{\omega }_{n,s},...,{\omega }_{N_{\mathrm{Assets}},N_T}\rangle .$

In this example, the quantum computing component which executes the parametrized quantum circuit PQC, prepares a state using two layers,

$\begin{array}{cc}❘\psi \rangle =U\left({\gamma }_2\right)U\left({\gamma }_1\right)❘{\psi }_0\rangle & \left(\mathrm{PQC}\right)\end{array}$

where the unitary operations U are freely parametrizable with a set of parameters γ={γ₁, γ₂}. The initial state of the circuit is given by the equal superposition of all possible states,

$❘{\psi }_0\rangle =\frac{1}{\sqrt{d_H}}\sum _x❘x\rangle ,$

where d_H=L^NAssetsNT is the dimensionality of the Hilbert space for N_Assets assets, N_T time slots and where the share of each asset can bet between zero and L−1, which are encoded by N_Assests·N_T qudits with local dimension L. In the current example, the unitary operation in each layer l=1, 2 has the form

$U\left({\gamma }_l\right)=\prod _{\alpha =1}^2{e}^{-i{\gamma }_{X;l,\alpha }X}{e}^{-i{\gamma }_{H;l,\alpha }{\mathscr{H}}_{\alpha }}$

which is a product over two different terms and each product term consists of a so-called mixing term ˜e^{−iγX;l,αX} and a phase operator ˜.

The mixing unitary has freely choosable parameters γ_X;l,α∈ and it is generated by the Hermitian mixing operator X. The mixing operators in this example is given by the sum over all x-angular momentum operators of each qudit, i.e.,

$X=\sum _{s=1}^{N_T}\sum _{n=1}^{N_{\mathrm{Assets}}}{L}_{x;n,s},\left(\mathrm{MIXEX}\right)$

where L_x;n,s is the usual x-angular momentum operator of the qudit encoding share of asset n at time slot s. The action of L_x;n,s is given by the standard action of angular momentum operators,

$L_{x;n,s}|{\omega }_{1,1},\dots ,{\omega }_{n,s},\dots ,{\omega }_{N_{\mathrm{Assets}},N_T}\rangle =\frac{1}{2}\sqrt{\left(L-{\omega }_{n,s}-2\right)\left({\omega }_{n,s}+2\right)}|{\omega }_{1,1},\dots ,{\omega }_{n,s},\dots ,{\omega }_{N_{\mathrm{Assets}},N_T}\rangle +\frac{1}{2}\sqrt{\left(L-{\omega }_{n,s}\right){\omega }_{n,s}}|{\omega }_{1,1},\dots ,{\omega }_{n,s},\dots ,{\omega }_{N_{\mathrm{Assets}},N_T}\rangle$

and as such generates transitions between and creates superpositions of basis states, i.e., different solutions.

Each of the two phase terms in each layer comes with one freely chooseable parameter, γ_H;l,α∈ and each term is generated by one cost Hamiltonians of (HAM1) and (HAM2).

Putting it all together, the quantum circuit prepares the following state,

$|\psi \left({\left\{{\gamma }_{X;l,1},{\gamma }_{H;l,1},{\gamma }_{X;l,2},{\gamma }_{H;l,2}\right\}}_{l=1}^2\right)\rangle =\prod _{l=1}^2{e}^{-i{\gamma }_{X;l,2}X}{e}^{-i{\gamma }_{H;l,2}{\mathscr{H}}_2}{e}^{-i{\gamma }_{X;l,1}X}{e}^{-i{\gamma }_{H;l,1}{\mathscr{H}}_1}❘{\psi }_0\rangle \left(\mathrm{PQCEX}\right)$

This parametrized quantum state depends on 8 parameters {γ_X;l,1, γ_H;l,1, γ_X;l,2, γ_H;l,2}_l=1² which need to be optimally determined by a classical optimization algorithm.

The proposed method and apparatus executes the following steps:

First, a classical optimization algorithm such as a stochastic gradient descent SGD is initiated on the classical computing component. The search parameters of this classical optimizer are given by the 8 parameters of the PQC, i.e., γ={γ_X;l,1, γ_H;l,1, γ_X;l,2, γ_H;l,2}_l=1²∈^D with D=8. Additionally, an archive is initiated on the classical computing component, where search parameters, candidate solutions and corresponding values of the objective functions will be stored.

Second, the quantum computing component is configured such that the measurement basis represents classical solutions as depicted in (ENC) and the quantum computing component produces the parametrized quantum state (PQCEX) with the corresponding gate operations (MIXEX) and (PQCEX) with (HAM1) and (HAM2). Then, an iteration counter is initialized with t=0 and the following steps are iterated:

• • i) Based on its internal state, the classical optimizer algorithm generates a set of search parameter values,

→γ^t={γ_X;l,1^t, γ_H;l,1^t, γ_X;l,2^t, γ_H;l,2^t}_l=1²

and the classical computing component sends these parameter values to the quantum computing component.

• • ii) The quantum computing component receives the parameter set

${\gamma }^t={\left\{{\gamma }_{X;l,1}^t,{\gamma }_{H;l,1}^t,{\gamma }_{X;l,2}^t,{\gamma }_{H;l,2}^t\right\}}_{l=1}^2$

and adjusts its unitary gate operators such that it produces the quantum state specified by (PQCEX), i.e.

$|\psi \left({\gamma }^t\right)\rangle =\prod _{l=1}^2{e}^{-i{\gamma }_{X;l,2}^{t}X}{e}^{-i{\gamma }_{H;l,2}^t{\mathscr{H}}_2}{e}^{-i{\gamma }_{X;l,1}^{t}X}{e}^{-i{\gamma }_{H;l,1}^t{\mathscr{H}}_1}❘{\psi }_0\rangle ,$

when it is executed.

Then the quantum apparatus executes the quantum circuit. The one resulting quantum state |ψ(γ^t) then encodes a whole population of N_S candidate solutions as superposition,

$|\psi \left({\gamma }^t\right)\rangle =\sum _{m=1}{\xi }_t\left(x_m\right)❘x_m\rangle .$

For this state a full state tomograph is done on the quantum device (see e.g. “Efficient quantum state tomography”, Nat Commun 1, 149 (2010) doi 10.1038/ncomms1147), which gives the list of all weights contributing to the state, T={|ξ_t(x_c)|² for c in 1, . . . d_H} and the set of N_S solutions _t at iteration t is extracted from this by selecting the ones which have the largest amplitudes, _t={x_m}_m=1^NS where x_m correspond to the largest N_S entries in the tomography list T.

• • iii) The classical computing component receives the population of solutions _t={x_m}_m=1^NS from the quantum computing component and stores it for further processing in the archive. Then, all objective function values from equations (OBJ1) and (OBJ2) are computed for each of the solutions in the population, which results in a set of objective value vectors _t={{right arrow over (C)}(x_m)∀x_m∈_t}.

From this set of objective value vectors, the set of non-dominated solutions _t={{right arrow over (C)}(x_i)}_i=1^NPF are selected. This represents the approximation of the Pareto front _t at iteration t and the corresponding set of solutions is an estimate of the Pareto set _t={x_i}_i=1^NPF. Form this set of nondominated solutions, the hypervolume indicator HV(_t) is calculated, which in turn is used to calculate the fitness value f_t associated with the parameters of the PQC,

$f_t\equiv f\left({\gamma }^t\right)=-HV\left({\mathrm{\mathbb{P}𝔽}}_t\right).$

The fitness value is then fed back into the classical optimization algorithm

←f(γ^t)

and the optimizer updates its internal state accordingly.

Additionally, the search parameters γ^t, the corresponding fitness value f(γ^t), the current estimate of the Pareto front _t, and the corresponding Pareto set of candidate solutions _t, are added to the archive,

→+{t, γ^t, f(γ^t), _t, _t}

• • iv) Thereafter, a convergence criterion is checked on the classical computation apparatus. This can be, for example, whether the maximally allowed number of cost function evaluations has been reached, or if the improvement of the fitness value f_t over the last couple of iterations was not enough.
    • a. If the convergence criterion is met the optimization is stopped and as the final result the final Pareto front and corresponding Pareto set are obtained by selecting all the non-dominated solutions from all approximations of the Pareto front of each iteration, , _Final=NonDominatedSolutions({_t′}_t′=1^t).
    • b. If the convergence criterion is not met, the process is repeated from step 1. above.

The final result is the Pareto set of solutions _Final={x_i}_i=1^NPF which approximate the Pareto front _Final={{right arrow over (C)}(x_i)}_i=1^NPF, which is then output and visualized. The Pareto front and the corresponding Pareto set form the basis for a financial service company to offer their customers various different portfolios as products. Each portfolio represents an optimal tradeoff solution between expected risk and expected return, and being able to offer solutions from all parts of the complete Pareto front as products to their customer give the financial service company a great benefit.

The proposed method and hybrid quantum-classical computer system utilizing the PQC is effective in producing a good approximation to the Pareto front since all possible solutions are typically present in the superposition of the initial state of the quantum apparatus. The phase operators in the PQC do assign different phases to each individual basis state based on the values of the scalarizing functions uses in the phase operators of Eq. (PQCEX). In that way, the classical optimizer can tune the parameters of the PQC to give the largest diversity and include mostly non-dominated solutions in the quantum state, which maximize hypervolume indicator.

Claims

1. A computer implemented method for obtaining a set of Pareto-optimal solutions to a given multi-objective optimization problem with M different objectives {right arrow over (C)}(x)={Cα}α=1M=(C1(x), C2(x), . . . , CM(x))T where x represents a classical solution of the multi-objective optimization problem, comprising performing an iterative procedure wherein each iteration t comprises the following steps: a. generating a set of parameters Yt on a classical computation component using the classical optimization algorithm and sending it to a quantum component,b. on the quantum component, generating a quantum state |ψγt=U(γt)|ψ0 by executing a parameterized quantum circuit (PQC) implementing a unitary quantum operation U(γ) with the set of parameters γt on initial quantum state |ψ0 where the quantum state |ψγt represents a set of NS different solutions,c. selecting a set of NS≥NPF basis states {|xi}i=1NS from the one quantum state |ψγt and sending the information of which basis states were selected to the classical computation component,d. on the classical computation component calculating a set of objective values consisting of the objective values for each solution in the set of classical solutions,e. generating a new estimate for the Pareto set of solutions and the corresponding Pareto front by the classical optimizer incorporating the set of solutions and the corresponding set of objectives values,f. updating the internal state of the classical optimization algorithm based on the current set of solutions and the Pareto set and their corresponding sets of objective values,g. checking a convergence criterion and if the convergence criterion is not satisfied, repeating the process from step a. onwards, until the convergence criterion is fulfilled,

2. The method according to claim 1, characterized in thatin addition to updating the internal state of the classical optimization algorithm, an archive is updated based on the current set of solutions and the Pareto set and their corresponding sets of objective values, in the step of generating the new estimate further incorporates values from previous iterations retrieved from an archive.

3. The method according to claim 1, characterized in thatthe parameterized quantum circuit PQC implements a unitary quantum operation U(γ) which is composed of NL layers with different choosable sets of parameters as

4. The method of claim 1, characterized in thatthe parameterized quantum circuit PQC implements a unitary quantum operation U(γ) which is composed of NL layers with different choosable sets of parameters as

5. The method according to claim 4, characterized in that the coefficients in these sums ak,α∈ represent reference vectors {right arrow over (a)}k which have to fulfill Σαak,α2={right arrow over (a)}kT·{right arrow over (a)}k=1 and ak,n≥0 for each k, and for each k, define a certain direction in the space of objective values and thus reflect one user preference and are determined in the beginning as part of the problem definition, αCα|xxe−iXk(γX;l,k)γX;l,kMXα(γX;l,k)=1, . . . , NP.

6. The method according to claim 1, characterized in thatthe parameterized quantum circuit PQC implements a unitary quantum operation U(γ) which is composed of NL layers with different choosable sets of parameters as

7. The method according to claim 1, characterized in thatthe basis-state selection is performed by sampling from the quantum state |ψγt which is realized by measuring repeatedly the quantum state, and adding for each measured state |xm the corresponding classical solution xm to the set of selected states, t until NS different states have been measured and added to the set of selected states, t={xi}i=1NS.

8. The method according to claim 1, characterized in thatthe basis-state selection is performed by a full state tomography from the quantum state |ψγt where the set of weights of each basis state is determined, T={|ξi|2}i=1dH, and the basis states corresponding to the subset of the NS largest weights is selected, which form the set of selected states, t={xi}i=1NS.

9. The method according to claim 3, characterized in thatthe mixing operator in each layer is given by a sum over all x-angular momentum operators or generalized Pauli X-operators of each physical qudit and where the coefficient of each individual angular momentum operator in each layer is given by one freely choosable parameter γX;k/α,l,i∈, where k/α enumerates the mixing operators, l denotes the layers and i enumerates all physical qudits.

10. The method according to claim 3, characterized in thatthe mixing operator in each layer is given by a sum over all x-angular momentum operators or generalized Pauli X-operators of each physical qudit and where there is only one global coefficient for the complete sum over all individual angular momentum or generalized Pauli X operators applied to each physical qudit in each layer which is given by the freely choosable parameter γX;k/α,l∈ where k/α enumerates the mixing operators, and l denotes the layers.

11. The method according to claim 4, characterized in thatthe coefficients ak,α of the reference vectors of the phase Hamiltonians, Pk=Σα=1Mak,αα, are determined and adapted in each iteration step t according to the state of the art procedure for reference vector adaptation such that in each iteration a different set of reference vectors is used which allows a speed-up of the optimization process and a diversification of the obtained Pareto solutions.

12. The method according to claim 1, characterized in thatthe problem is a multi-objective bi-directional electric vehicle charging problem where the target is to find a set of Pareto optimal configurations of discrete integer variables, which encode a full charging and service schedule for a fleet of electric vehicles and where multiple objectives, which are at least two of electricity cost for fulfilling all charging demands, customer satisfaction, maximizing utilization of renewable energies, minimizing electricity drawn from the grid, and minimizing the peak electricity power consumption, are considered for achieving an approximation of the Pareto front, is obtained.

13. The method according to claim 1, characterized in thatthe method is applied to obtain the efficient frontier of possible portfolio compositions as formulated by the Markowitz model of Modern portfolio theory by finding the compositions of a portfolio over time as specified by integer variables which specify the investment amount into asses n at time s and where the efficient frontier is spanned by those portfolio compositions which have optimal tradeoffs between maximizing return and minimizing risk and thus pose a two-objective optimization problem.

14. A hybrid quantum-classical computer system comprising a quantum component configured to realize a parametrized quantum circuit PQC by implementing a unitary quantum operation U(γ) which is parametrized by freely choosable parameters γ and which generates a N-qubit or qudit quantum state |ψγ by applying the unitary operator to a previously defined N-qubit or qudit initial state |ψ0, i.e. |ψγ=U(γ)|ψ0, which can be decomposed into the measurement basis of the quantum device, |ψγ=Σc=1dHξγ,c|xc with dH the dimensionality of the Hilbert space, and each measurement basis state |xc encodes one classical solution xc,and a classical computing component which is configured to run a classical optimization algorithm by providing new parameters values γ to the quantum computing component, initiating the execution of the PQC and the preparation of the state |ψγ and then receiving a set of classical solutions from the quantum module as a result and calculating the all cost functions for each solutions,where both modules together are used to solve classical multi-objective optimization problem by iteration as defined in the method according to claim 1.