Controller for early termination in mixed integer optimal control optimization

Abstract

A hybrid integer optimal control optimization problem is solved by using a branch and bound (B&B) optimization that searches for a global optimal solution within a search space to control a system. The B&B optimization iteratively partitions the search space into a region-nested tree and prunes at least one region from the region-nested tree before finding a local optimal solution for each region when a dual objective value of a projection of a suboptimal dual solution estimate for each region in a dual feasible space is greater than an upper bound of a global optimal solution maintained by the B&B optimization or less than a lower bound of the global optimal solution.

Description

Technical Field

[0001] This disclosure generally relates to optimization-based control, and more specifically to methods and apparatus for model predictive control of systems having continuous and discrete operational elements. Background Technology

[0002] Optimization-based decision-making, planning, and control techniques, such as Model Predictive Control (MPC), allow for model-based design frameworks in which system dynamics, requirements, and constraints can be directly considered. This framework has been extended to include hybrid systems with both continuous and discrete decision variables, providing powerful techniques for modeling a wide range of problems, such as dynamic systems with mode switching or quantized actuation, logical rule problems, time-logic rules, or obstacle avoidance constraints. However, the resulting optimization problems are highly non-convex, making them difficult to solve in practice due to the inclusion of variables that only take integer values. When linear system dynamics and linear inequality constraints are combined with linear or linear quadratic objectives, the resulting optimal control problem (OCP) can be formulated as Mixed Integer Linear Programming (MILP) or Mixed Integer Quadratic Programming (MIQP). More general convex inequality constraints, such as quadratic inequality constraints, can be included, resulting in Mixed Integer Quadratic Constraint Quadratic Programming (MIQCQP), or second-order cone constraints, resulting in Mixed Integer Second-Order Cone Programming (MISOCP).

[0003] Decision making, planning, or control of hybrid systems aims to solve a mixed-integer programming problem (MIP) at each sample time. This is a challenging task because MIPs are typically NP-hard, and several approaches to solving this sequence of MIPs have been explored in the literature. These approaches can be categorized into heuristic techniques seeking efficient suboptimal solutions to the problem and optimization algorithms attempting to solve the MIP to optimality. Examples of the former techniques can be based on rounding and pumping schemes using approximate optimization algorithms, or on approximate dynamic programming. A drawback of fast heuristics is that they typically cannot guarantee finding an optimal or even integer-feasible but suboptimal solution.

[0004] Most mixed-integer optimization algorithms are based on variations of the branch and bound (B&B) technique to bring MIP solutions to optimality. Variations of the branch and bound strategy have been combined with various methods for solving relaxed convex subproblems, such as dual active set solvers, interior-point algorithms, dual projected gradient methods, non-negative least squares solvers, and alternating direction multipliers (ADMM). One advantage of using dual active set solvers is that the dual feasible starting point of the subproblem can be computed from the dual solution of the parent problem. The latter allows the solver of the relaxed problem to terminate early when the biobjective becomes greater than the current upper bound, thus reducing the computational cost of branch and bound optimization methods.

[0005] Interior-point methods (IPMs) can be used to solve more general convex optimization problems, and they typically exhibit better worst-case computational performance compared to active-set optimization algorithms. Furthermore, recent advances in warm-up strategies for IPMs have spurred their use in branch-and-bound optimization of mixed-integer programming. However, unlike dual active-set methods, IPMs typically do not provide a dual feasible solution estimate in each iteration, which can be used in convex solvers that terminate branch-and-bound strategies early. Summary of the Invention

[0006] This disclosure proposes a method for early termination before computing the optimal solution to a convex relaxation, based on the detection that a solution to a convex relaxation might be discarded due to infeasibility or based on the current global bound, in order to reduce the computational workload of the branch-and-bound method. This disclosure describes a custom projection step to compute a dual feasible solution estimate given the primal dual solution estimate in the current iteration of the interior-point method. This allows for evaluation of the dual objective used to early terminate the convex solver in the branch-and-bound method. The same early termination strategy can be used for primal infeasibility detection.

[0007] The implementation method is based on the branch and bound (B&B) optimization algorithm for solving mixed-integer programming problems. The B&B optimization algorithm for mixed-integer programming requires solving multiple convex relaxations to calculate the lower and upper bounds of the optimal solution for the mixed-integer programming problem. When a convex relaxation is originally infeasible or when the objective value exceeds the current global upper bound, multiple solutions with convex relaxations can usually be discarded. Interior-point methods are generally applicable to a wide range of optimization problems and can be used very reliably to solve convex relaxations in branch and bound methods for mixed-integer programming, even without a good initial estimate of the optimal solution.

[0008] Therefore, some implementations aim to use B&B optimization algorithms to solve mixed-integer programming problems for decision making, planning, or control of mixed systems involving continuous and discrete decision variables. Additionally or alternatively, another objective of some implementations is to provide a control technique using the interior-point method (IPM) of convex relaxation in branch-and-bound methods for solving mixed-integer programming problems.

[0009] Some implementations are based on the understanding that the advantage of B&B optimization lies in its ability to effectively reduce the search space by eliminating spaces that do not contain optimal results. Specifically, the goal of the B&B algorithm is to find a value x within a set S of acceptable or candidate solutions that maximizes or minimizes the value of a real-valued function f(x), called the objective function. This set S is called the search space or feasible region. Furthermore, S... I The set of candidate solutions used to represent instance I. Instance representation includes three operations such as branching (I), bounding (I), and solving (I). The branching operation (I) produces each representative S. I Two or more instances of a subset. Typically, subsets are disjoint to prevent the algorithm from visiting the same candidate solution twice, but this is not required. However, S I The optimal solution must be contained in at least one subset. Operational bounding (I) computes a lower bound (and / or upper bound) on the value of any candidate solution in the space represented by I, i.e., for S I For all x, the bound (I) ≤ f(x). The operation solves (I) to determine whether I represents a single candidate solution. Alternatively, if not, the operation can choose to return S. I Some of the feasible solutions.

[0010] Using these operations, the B&B algorithm performs a top-down recursive search through a tree of instances formed by the branching operations. Once an instance I is visited, it checks if the bound (I) is greater than the lower bound and / or greater than the global upper bound of some other instance that has already been visited. If the bound (I) is greater than the global upper bound, I can be safely discarded from the search and recursive steps. This pruning step is typically implemented by maintaining a global variable that records the minimum lower bound and maximum upper bound seen among all instances checked so far.

[0011] Typically, to evaluate the lower bound of instance I representing a partition of the search space, it is necessary to determine a solution for a convex relaxation of instance I. In many practical applications of B&B algorithms for solving non-convex mixed integer problems constrained by equality and inequality, the calculation of the lower bound of instance I is performed using a convex relaxation of the optimization problem on the search space corresponding to instance I. Since many convex relaxations often need to be solved as part of the B&B algorithm to find the optimal solution to a single mixed integer problem, it is important to efficiently solve for convex relaxations and to reduce the computational burden of the B&B algorithm.

[0012] Some implementations are based on the understanding that, in optimization theory, duality, or the duality principle, is a principle that allows us to view optimization problems from either two perspectives (formulating the primal or dual problem). The solution to the dual problem provides a lower bound on the solution to the primal (minimizing) problem. Typically, the optimal values ​​of the primal and dual problems are not necessarily equal. Their difference is called the duality gap. However, for many convex optimization problems, the duality gap is zero under constraints.

[0013] Some implementations are based on the understanding that the solution to the dual problem provides a lower bound on the solution to the primal (minimizing) problem. If the dual objective value of a feasible but suboptimal solution to the dual problem is greater than an upper bound previously determined in the B&B algorithm used to solve mixed-integer programming, then the optimal solution to the dual problem is also greater than the upper bound, and any feasible solution to the primal problem is also greater than the upper bound. Therefore, in the context of the B&B method, if the dual objective value of at least one feasible solution to the dual problem is greater than the upper bound, it is not necessary to find the optimal solution to the dual problem and / or the primal problem. In other words, if the dual objective value of any dual feasible solution to instance I is greater than the upper bound maintained by the B&B method, then the optimal solution to instance I cannot be lower than the upper bound and should be discarded. In this way, the computational burden of the B&B method can be reduced.

[0014] However, some implementations are based on the following understanding: unfortunately, these principles for early termination of the convex relaxation solution process based on dual feasible solution estimation cannot be directly applied when evaluating solutions to convex relaxations using the primordial-dual interior-point method. This is because the primordial-dual interior-point method is designed to find both feasible and optimal solutions simultaneously, meaning that the primordial or dual feasibility of intermediate solutions cannot be guaranteed until the optimal solution is found. Early termination based on the dual objective requires dual feasible solution estimation.

[0015] Some implementations are based on the understanding that primal dual interior-point optimization is an iterative method. The feasibility of intermediate and suboptimal dual solutions can be tested against each other, and the suboptimal dual solution is compared to the upper bound only if it is feasible. However, due to the principles of interior-point optimization, the feasibility of suboptimal dual solutions is more likely to lie in a space closer to optimality, thus the computational cost of finding the estimate of the dual feasible solution in standard interior-point optimization methods can be substantial. Therefore, testing feasibility and testing the dual objective value of the estimate of the dual feasible solution against the upper bound to eliminate suboptimal branches in the B&B algorithm is more likely to be close to the end of the interior-point optimization method. However, it is desirable to eliminate suboptimal branches as quickly as possible to reduce the computational burden of the B&B method.

[0016] Some implementations are based on the understanding that instead of checking the feasibility of the suboptimal dual solution, it is more efficient to project the suboptimal dual solution into its feasible space and compare the projected dual objective value with the upper bound. In this way, even infeasible solutions can be used to check the optimality of the current evaluation instance I without finding its optimal dual solution and / or primal solution.

[0017] Some embodiments of this disclosure are based on the interior point method (IPM) to solve convex relaxation of inequality constraints in branch and bound algorithms, and the implementation of the proposed method includes a customized projection step based on the solution of the equality constraint optimization problem, which can be used to compute a dual feasible solution estimate for early termination of the solver, starting from the primal dual solution estimate of the current iteration of the IPM.

[0018] In interior-point optimization methods, for inequality constraints that are inactive at the optimal solution (i.e., when the left and right sides of such inequality constraints are not equal at the solution), the Lagrange multipliers become increasingly closer to zero while the slack variables remain non-zero. Conversely, for inequality constraints that are active at the optimal solution (i.e., when the left and right sides of such inequality constraints are equal at the solution), the slack variables become increasingly closer to zero while the Lagrange multipliers remain non-zero. Regarding the ratio between the slack variables and the Lagrange multipliers, in subsequent iterations of IPM, this ratio increases or decreases, respectively, for each inequality constraint that is inactive or active at the optimal solution with convex relaxation.

[0019] The embodiments of this disclosure are based on the understanding that solving the inequality-constrained optimization problem in the projection step to ensure the feasibility of both equality and inequality constraints with respect to the dual problem of convex relaxation is computationally expensive, and that solving the equality-constrained optimization problem is preferable to implementing the projection step. Some embodiments are based on the understanding that, given the primal-dual solution estimate in the current iteration of the IPM, the value of the ratio between the relaxation variables and the Lagrange multiplier values ​​can be used in the weighted matrix of the Lagrange multipliers in the objective of the equality-constrained optimization problem solved in the projection step that computes the dual feasible solution estimate. The solution to the latter equality-constrained optimization problem offers several important advantages compared to the standard minimum norm projection of dual feasibility constraints.

[0020] First, based on the values ​​of the slack variables The value of the Lagrange multiplier for each inequality constraint The ratio between Using the objective term in the proposed equality-constrained optimization problem This achieves the projection step, indirectly forcing a positive constraint μ≥0. More specifically, at the current value... When it approaches zero, the target item Penalty step size Δμ i Keep it small to get a ratio The large value of this makes it possible for it to satisfy the positive constraint of the Lagrange multiplier value after the projection step.

[0021] Secondly, in some embodiments of this disclosure, the optimal solution to the equality-constrained optimization problem in the proposed projection step is computed by solving a structured linear system corresponding to the first-order necessary optimality conditions of convex relaxation (i.e., the Karush-Kuhn-Tucker (KKT) conditions). In some embodiments, the linear system in the projection step can be solved using the same matrix factorization process used to solve the linear system in each iteration of IPM. Additionally, some embodiments reuse the matrix factorization from the projection step to solve the linear system in subsequent IPM iterations. Finally, unlike standard minimum norm projection, the equality-constrained optimization problem in the proposed projection strategy aims to preserve any progress toward IPM convergence, i.e., it preserves the original feasibility and optimality of the primal dual solution estimate for IPM.

[0022] If the dual objective value is greater than a certain threshold, given the original dual solution estimate at the current iteration of the IPM, the implementation of this disclosure uses the proposed projection step to compute the dual feasible solution estimate. More specifically, if the dual objective value exceeds the current global upper bound in the branch-and-bound optimization method, then after the projection step, for the dual feasible solution estimate, the dual objective may exceed the upper bound, causing the solution of the convex relaxation IPM to terminate prematurely and the computational workload to be reduced accordingly.

[0023] The embodiments of this disclosure are based on the understanding that the objective value of the dual problem is less than or equal to the objective value of the original formula for convex relaxation, such that the dual objective value exceeds the current upper bound in the branch and bound method only when the optimal objective value of convex relaxation exceeds the current upper bound, that is, only when the solution of convex relaxation can be discarded.

[0024] When the primal problem corresponding to the convex relaxation is infeasible, the dual objective is unbounded. The implementation of this disclosure is based on the understanding that for subsequent iterations of the IPM applied to the primally infeasible convex relaxation, the dual objective reaches an infinite value, such that the proposed projection step can be used to compute an estimate of the dual feasible solution, and then prematurely terminates whenever the value of the dual objective exceeds the current upper bound in the branch and bound method.

[0025] Some implementations involve branch-and-bound optimization methods to solve mixed-integer optimal control problems (MI-OCP), which involve both continuous and discrete variables in the model, describing the system, performance indices, and system requirements, respectively, and the objectives and constraints. Some implementations include continuous or discrete state and control input variables; that is, they can only assume discrete sets of values. This disclosure aims to solve the block-sparse mixed-integer optimal control problem to achieve model predictive control of hybrid systems, enabling the use of problem sparsity in our proposed early termination strategy to efficiently solve the relaxation of the convex optimization problem in the branch-and-bound method.

[0026] In some embodiments of this disclosure, the projection step may use Cholesky factorization based on a combination of forward and backward substitutions to solve an equality-constrained optimization problem that provides updated Lagrange multiplier values ​​to achieve duality. In some embodiments, block tridiagonal Cholesky factorization may be used in conjunction with block-structured forward and backward substitutions. In some embodiments, the (block-structured) Cholesky factorization from the projection step may be reused to solve the linear system in subsequent IPM iterations.

[0027] Some implementations are based on the understanding that the computational cost of evaluating the dual objective once is significantly less than the computational cost of a projection step. Therefore, in embodiments of this disclosure, a projection step is performed if and only if the evaluation value of the dual objective estimated by the dual infeasible solution is greater than the current upper bound in the branch and bound optimization method.

[0028] Therefore, an embodiment discloses a predictive controller for a control system, the predictive controller comprising: at least one processor; and a memory having instructions stored thereon, the instructions, when executed by the at least one processor, causing the predictive controller to: receive a feedback signal including a measurement result of the system's state; solve a mixed-integer optimal control optimization problem using branch and bound (B&B) optimization that searches for a global optimum within a search space to generate a control signal, wherein the B&B optimization iteratively partitions the search space into a nested tree of regions, and removes at least one region from the nested tree before finding a local optimum for each region when the dual objective value of the projection of the suboptimal dual solution estimate of each region into the dual feasible space is greater than an upper bound of the global optimum maintained by the B&B optimization or less than a lower bound of the global optimum maintained by the B&B optimization; and control the system based on the control signal to change the state of the system.

[0029] According to the example implementation, the region nesting tree is formed by different convex relaxations of integer variables in a mixed integer optimal control optimization problem, and wherein a local optimum is searched in each region of the region nesting tree over multiple iterations by an interior-point method, such that a suboptimal dual solution estimate for each region is generated by intermediate iterations of the interior-point method.

[0030] According to the example implementation, the local optimum is based on an early termination process, which allows the interior-point method to terminate its iteration process before finding a local optimum in at least one region if: after one or more iterations of the interior-point method, a dual feasible solution estimate is computed, and the dual objective value is greater than the upper bound of the global optimum maintained by B&B optimization; or, after one or more iterations of the interior-point method, a projection step is performed to compute a dual feasible solution estimate, and the corresponding dual objective value is greater than the upper bound of the global optimum maintained by B&B optimization; or, after one or more iterations, the interior-point method detects that no feasible solution exists for at least one region.

[0031] According to the example implementation, the processor is configured to: when the dual objective value estimated by the suboptimal dual solution is greater than the upper limit, perform a projection of the suboptimal dual solution estimate for early termination of the process; and otherwise continue the iterative process from the suboptimal dual solution estimate to find the optimal solution of convex relaxation.

[0032] According to the example implementation, the projection step computes an update to the estimate of the suboptimal dual solution by solving an equality-constrained optimization problem that enforces dual feasibility constraints and minimizes the norm of the update to the Lagrange multiplier variables.

[0033] According to the example implementation, the local optimal solution of at least one region is used to initialize the iterative solution process in one or more other regions within the nested region tree for B&B optimization.

[0034] According to the example implementation, the projection step computes a suboptimal dual solution estimate from the local optimum to perform an early termination process for one or more other regions within the nested region tree used for B&B optimization.

[0035] According to the example implementation, the projection step imposes one or more additional optimality conditions based on the positive definite weighting matrix in the objective function and minimizes the weighted norm of the updated optimization variables.

[0036] According to the example implementation, the projection step calculates an update to the estimate of the suboptimal dual solution by solving the block-structured linear system, the solution of which forms the optimal solution to the equality-constrained optimization problem.

[0037] According to the example implementation, the positive definite weighting matrix in the objective of the projection step corresponds to a weighted Hessian matrix as an update of the original optimization variables, a diagonal matrix with small positive weights for updating the Lagrange multipliers corresponding to the equality constraints in the convex relaxation, and a positive weighting value corresponding to the ratio between the relaxation variables and the Lagrange multiplier values ​​corresponding to each inequality constraint after one or more iterations in the convex solution process.

[0038] According to the example implementation, the block-structured linear system is solved by solving the simplified linear system to compute the update of the optimization variables, and then solving the diagonal linear system to compute the update of the Lagrange multiplier variables.

[0039] According to the example implementation, given a Cholesky factorization from a prior projection step or a prior iteration from a convex solution process, the Cholesky factorization is computed or updated to compute a solution to a simplified linear system with a positive definite KKT matrix.

[0040] According to the example implementation, the KKT matrix of the simplified linear system has a block tridiagonal sparse structure due to the block-structured sparsity of the mixed-integer optimal control optimization problem, and the block tridiagonal Cholesky factor decomposition is computed or updated to compute the solution of the simplified linear system with a positive definite block tridiagonal KKT matrix.

[0041] According to the example implementation, the predictive controller is implemented using mixed-integer model predictive control (MI-MPC), wherein MI-MPC calculates the control signal based on the current state of the system and the control command, and wherein MI-MPC calculates the control solution by solving a constrained mixed-integer optimization problem at each control time step, the control solution comprising a sequence of future optimal discrete and continuous control inputs over the system's predictive time history.

[0042] According to the example implementation, the system is a vehicle, and the predictive controller determines the inputs to the vehicle based on a mixed-integer control solution, wherein the inputs to the vehicle include one or a combination of the vehicle's acceleration, engine torque, braking torque, and steering angle, and discrete optimization variables are used to model one or a combination of discrete control decisions, switching, gear shifting, and obstacle avoidance constraints in system dynamics.

[0043] According to the example implementation, the system is a spacecraft, and the predictive controller determines the input to the spacecraft based on a mixed integer control solution, wherein the input to the spacecraft drives one or a combination of thrusters and momentum exchange devices, and discrete optimization variables are used to model one or a combination of discrete control decisions, switching in system dynamics, integer values ​​of thruster commands, and obstacle avoidance constraints.

[0044] According to one example implementation, the system is a vapor compression system, and the predictive controller determines the inputs to the vapor compression system based on a mixed integer control solution, wherein the inputs to the vapor compression system include one or a combination of indoor unit fan speed, outdoor unit fan speed, compressor rotational speed, expansion valve position, and flow reversal valve position, and discrete optimization variables are used to model one or a combination of discrete control decisions, switching in system dynamics, and integer values ​​of commands sent to valves and / or fans.

[0045] Accordingly, another embodiment discloses a method for controlling a control system, wherein the method uses at least one processor; and a memory having instructions stored thereon, the instructions implementing the steps of the method when executed by the at least one processor, the steps including: receiving a feedback signal including a measurement result of the state of the system. The method further includes solving a mixed-integer optimal control optimization problem using branch-and-bound (B&B) optimization that searches for a global optimum within a search space to generate a control signal, wherein the B&B optimization iteratively partitions the search space into a nested tree of regions, and removes at least one region from the nested tree before finding a local optimum for each region if the dual objective value of the projection of the suboptimal dual solution estimate of each region into the dual feasible space is greater than an upper bound or less than a lower bound of the global optimum maintained by the B&B optimization. The method also includes controlling the system based on the control signal to change the state of the system.

[0046] According to an exemplary embodiment, a non-transitory computer-readable storage medium includes a processor-executable program for performing a method comprising: receiving a feedback signal including a measurement result of the state of a system; solving a mixed-integer optimal control optimization problem using branch and bound (B&B) optimization that searches for a global optimum within a search space to generate a control signal, wherein the B&B optimization iteratively partitions the search space into a nested tree of regions, and prunes at least one region from the nested tree before finding a local optimum for each region when the dual objective value of the projection of the suboptimal dual solution estimate of each region into the dual feasible space is greater than an upper bound or less than a lower bound of the global optimum maintained by the B&B optimization; and controlling the system based on the control signal to change the state of the system. Attached Figure Description

[0047] [ Figure 1A ]

[0048] Figure 1A Block diagrams of predictive controllers and feedback systems according to some implementation methods are illustrated.

[0049] [ Figure 1B ]

[0050] Figure 1B Block diagrams of predictive controllers and feedback systems according to some implementation methods are illustrated.

[0051] [ Figure 2A ]

[0052] Figure 2A A block diagram illustrating a method for a predictive controller, according to some implementations, for mixed-integer model predictive control (MI-MPC) to calculate control signals based on the current state of a given system and control commands is shown.

[0053] [ Figure 2B ]

[0054] Figure 2B A block diagram illustrating the MI-MPC method for solving optimal control structured mixed integer quadratic programming (MI-QP) according to some implementations is shown.

[0055] [ Figure 3A ]

[0056] Figure 3A A schematic diagram illustrating an example of a binary control variable search tree representing a nested tree of search regions for feasible control solutions of integers, according to some implementation methods.

[0057] [ Figure 3B ]

[0058] Figure 3B A block diagram is shown, according to some implementations, of a branch-and-bound mixed integer optimization algorithm for searching for integer feasible optimal control solutions based on a nested tree of search regions and corresponding lower / upper bounds.

[0059] [ Figure 4A ]

[0060] Figure 4A The original and dual formulas of convex relaxation, which are solved as part of the branch-bound optimization algorithm in the MI-MPC controller according to some implementations, are illustrated.

[0061] [ Figure 4B ]

[0062] Figure 4B The connection between the primal and dual formulas for each convex relaxation in a branch-and-bound optimization algorithm according to some implementations is illustrated.

[0063] [ Figure 4C ]

[0064] Figure 4C This illustrates the solution process for early termination of convex relaxation using suboptimal dual solution estimation within a branch-and-bound optimization method, based on dual feasible solution estimation or projection onto dual feasible solution estimation, according to some implementations.

[0065] [ Figure 5A ]

[0066] Figure 5A A block diagram illustrating an early termination process for solving convex relaxations within a branch-and-bound optimization method, based on dual objective values ​​estimated from dual feasible solutions, is shown according to some implementations.

[0067] [ Figure 5B ]

[0068] Figure 5B A block diagram illustrating an early termination process for solving convex relaxations within a branch-and-bound optimization method, based on a projection step that estimates dual feasible solutions, is presented according to some implementations.

[0069] [ Figure 6A ]

[0070] Figure 6A An example is given of a block sparse interior-point optimization algorithm for solving constrained optimal control structured quadratic programming (OCP-QP) by iteratively solving a relaxed system with necessary optimality conditions, according to some implementations.

[0071] [ Figure 6B ]

[0072] Figure 6B A block diagram illustrating the initialization steps and iterative process of an interior-point optimization algorithm for solving constrained optimal control problems in the implementation of predictive control systems, according to some implementations, is provided.

[0073] [ Figure 6C ]

[0074] Figure 6C A block diagram illustrating the computation of a linearized KKT system with a Newton-type search direction during the iterative process of solving an interior point optimization algorithm for constrained OCP-QP, according to some implementations, is shown.

[0075] [ Figure 6D ]

[0076] Figure 6D A block diagram illustrating a sequence of decompositions of a linearized KKT system into two subsystems in an interior-point optimization algorithm for solving constrained OCP-QP, according to some implementations, is shown.

[0077] [ Figure 7A ]

[0078] Figure 7A A block diagram illustrating the process of calculating the projection step within a branch-and-bound optimization algorithm, which can be used to prematurely terminate the estimation of the dual feasible solution of a convex solver, according to some implementations.

[0079] [ Figure 7B]

[0080] Figure 7B A block diagram illustrates the process of calculating the projection steps of the dual feasible solution estimate by solving an equality-constrained optimization problem to enforce affine equality constraints in a convex relaxed dual QP formula, according to some implementations.

[0081] [ Figure 7C ]

[0082] Figure 7C A block diagram illustrates the process of calculating the projection step of the dual feasible solution estimate by solving an equality-constrained linear quadratic optimization problem in the objective of the equality-constrained QP with alternative scaling of the optimization variables to effectively enforce the positivity of the updated Lagrange multiplier values.

[0083] [ Figure 7D ]

[0084] Figure 7D Pseudocode is illustrated for solving the early termination process of convex relaxation in a branch-and-bound optimization algorithm in a predictive controller, according to some implementations, for evaluation based on the dual objective and projection steps based on the estimation of dual feasible solutions.

[0085] [ Figure 8 ]

[0086] Figure 8 An example is given of the equivalent of the equation-constrained QP solution and the solution of the simplified linear system in the computational projection step according to some implementations, and then the equations of the approximately linearized KKT system with dual variable updates are independently computed.

[0087] [ Figure 9A ]

[0088] Figure 9A A block diagram illustrates a process for solving a simplified approximate linearized KKT system based on a forward Cholesky factorization of a positive definite approximate KKT matrix, according to some implementations, to compute the Newton-type search direction of the original variables in order to compute the projection step of the Lagrange multipliers.

[0089] [ Figure 9B ]

[0090] Figure 9B A block diagram illustrates a process for solving a simplified approximate linearized KKT system by inverse Cholesky factorization based on a positive definite approximate KKT matrix, according to some implementations, to compute the Newton-type search direction of the original variables.

[0091] [ Figure 9C ]

[0092] Figure 9C This is an example of a block-structured sparse pattern diagram of an approximate KKT matrix in a linearized system of Newton-type search directions used to compute the original and dual optimization variables in the projection step, according to some implementations.

[0093] [ Figure 9D ]

[0094] Figure 9D This illustrates a block structure sparse pattern diagram of a positive (semi-)definite KKT matrix in a simplified approximate linearized KKT system with a Newton-type search direction for the original variables in the projection step, according to some implementations.

[0095] [ Figure 9E ]

[0096] Figure 9E A block diagram illustrates a process for solving a simplified approximate linearized KKT system based on a block triangular forward Cholesky factorization of a block tridiagonal positive definite approximation KKT matrix, according to some implementations, to compute the Newton-type search direction of the original variables in order to compute the projection step of the Lagrange multipliers.

[0097] [ Figure 9F ]

[0098] Figure 9F A block diagram illustrates the process of solving a simplified approximate linearized KKT system by block tridiagonal inverse Cholesky factorization based on a block tridiagonal positive definite approximation KKT matrix to compute the Newton-type search direction of the original variables, according to some implementations.

[0099] [ Figure 10A ]

[0100] Figure 10A A block diagram illustrates a warm-start initialization process for calculating the solution estimate of a convex solution process based on an approximate or optimal solution of a convex relaxation in the parent node, according to some implementations.

[0101] [ Figure 10B ]

[0102] Figure 10B A block diagram illustrating a smooth relaxation-based warm-start initialization process for solving one or more convex relaxation interior point optimization algorithms in a predictive control system, according to some embodiments, is shown.

[0103] [ Figure 11 ]

[0104] Figure 11A block diagram illustrating an early termination process for solving convex relaxation within a branch and bound algorithm, according to some implementations, is shown. This process is used when the optimal objective value is greater than the current upper limit or when no feasible solution exists for convex relaxation.

[0105] [ Figure 12A ]

[0106] Figure 12A A schematic diagram of a vehicle is shown, illustrating a predictive controller that employs principles based on some implementation methods.

[0107] [ Figure 12B ]

[0108] Figure 12B A schematic diagram illustrating the interaction between a vehicle's predictive controller and other controllers according to some implementations is shown.

[0109] [ Figure 12C ]

[0110] Figure 12C A schematic diagram illustrating a path and / or motion planning method for a controlled vehicle, employing principles from some implementation methods, is shown.

[0111] [ Figure 12D ]

[0112] Figure 12D Exemplary traffic scenarios are illustrated based on single-vehicle or multi-vehicle decision-making modules according to some implementation methods.

[0113] [ Figure 13A ]

[0114] Figure 13A It is a schematic diagram of the formulation of the spacecraft mixed integer predictive control problem based on the principles of some implementation methods.

[0115] [ Figure 13B ]

[0116] Figure 13B It is a schematic diagram of the formulation of the spacecraft mixed integer predictive control problem based on the principles of some implementation methods.

[0117] [ Figure 14A ]

[0118] Figure 14A A schematic diagram of a controller-controlled vapor compression system according to some embodiments is shown.

[0119] [ Figure 14B ]

[0120] Figure 14B Examples of signal, sensor, and controller configurations used in a VCS according to some implementations are illustrated.

[0121] [ Figure 15 ]

[0122] Figure 15 A method for controlling a system according to an example implementation is illustrated. Detailed Implementation

[0123] Some embodiments of this disclosure provide systems and methods for operating a control system or a system using a predictive controller. An example of a predictive controller is model predictive control (MPC), which determines the control input based on a model of the controlled system.

[0124] Figure 1A A block diagram illustrating a predictive controller 110 and a feedback system 120 according to some embodiments is shown. Figure 1A An example feedback system (or system) 120, according to some embodiments, is connected to a predictive controller 110 (or controller) via a state estimator 130. In some implementations, the predictive controller 110 is an MPC controller programmed according to a dynamic model 102 (or system model) of the system 120. The system model 102 may be a set of equations representing the system 120 as a function of current and previous inputs 111 and previous outputs 103, showing how the state and output 103 change over time. The system model 102 may include constraints 104 representing physical and operational limitations of the system 120. During operation, the controller 110 receives a command 101 representing the desired behavior of the system 120. The command may be, for example, a motion command. In response to receiving the command 101, the controller 110 generates a control signal 111 to serve as an input to the system 120. In response to the input, the system updates the output 103 of the system 120. Based on measurements of the output 103 of the system 120, the estimator 130 updates the estimated state 121 of the system 120. The estimated state 121 of system 120 provides state feedback to predictive controller 110. Therefore, the predictive controller initially receives the feedback signal 121 from system 120 via estimator 130, wherein the feedback signal 121 includes the measurement result of the state of system 120.

[0125] As referred to herein, system 120 can be any machine or device controlled by certain manipulating input signals (e.g., control signal 111 (input)) that may be associated with physical quantities such as voltage, pressure, force, torque, and return certain controlled output signals 103 (output) that may be associated with physical quantities such as current, flow, velocity, position, etc., indicating a transition of the system state from a previous state to the current state. The output value is partly related to the system's previous output value and partly related to both the previous and current input values. The dependence on previous inputs and outputs is encoded in the system state. The operation of the system (e.g., the movement of system components) can include a sequence of output values ​​generated by the system after certain input values ​​are applied.

[0126] System model 102 may include a set of mathematical equations describing how the system output changes over time as a function of the current input, previous inputs, and previous outputs. The state of system 120 is a collection of any information that typically changes over time, such as the current input and appropriate subsets of previous inputs and outputs, which, together with the system model and future inputs, can uniquely define the future motion of the system.

[0127] System 120 may be subject to physical limitations and specification constraints 104, which limit the range of operations that system 120's outputs, inputs, and possible states are allowed to operate.

[0128] The controller 110 may be implemented in hardware or as a software program executing in a processor (e.g., a microprocessor), which receives the estimated state 121 of the system 120 and the desired motion command 101 at fixed or variable control cycle sampling intervals and uses this information to determine inputs for the operating system 120, such as control signals 111. The controller 110 also uses branch and bound (B&B) optimization, which searches for a global optimum within the search space, to solve a mixed-integer optimal control optimization problem to generate the control signals 111. B&B optimization iteratively partitions the search space into a nested tree of regions, and removes at least one region from the nested tree before finding a local optimum for each region if the dual objective value of the projection of the suboptimal dual solution estimate in the dual feasible space of each region is greater than the upper bound of the global optimum maintained by B&B optimization or less than the lower bound of the global optimum maintained by B&B optimization. Furthermore, the nested tree of regions is formed by different convex relaxations of the integer variables of the mixed-integer optimal control optimization problem. The controller 110 also controls the system 120 based on the control signals 111 to change the state of the system 120.

[0129] Furthermore, by searching for local optima within each region of the nested region tree in multiple iterations using the interior-point method, suboptimal dual solution estimates for each region are generated through intermediate iterations of the interior-point method. Local optima are based on an early termination process, which allows the interior-point method to terminate its iteration process before finding a local optimum in at least one region if:

[0130] 1. A dual feasible solution estimate is calculated after one or more iterations of the interior-point method, and the dual objective value is greater than the upper bound of the global optimum maintained by B&B optimization; or

[0131] 2. Perform a projection step after one or more iterations of the interior-point method to compute an estimate of the dual feasible solution, and the corresponding dual objective value is greater than the upper bound of the global optimum maintained by B&B optimization; or

[0132] 3. The interior point method detects after one or more iterations that there is no feasible solution for at least one region.

[0133] Furthermore, the local optimal solution for at least one region of the region nesting tree is used to initialize the iterative process in one or more other regions within the region nesting tree used for B&B optimization.

[0134] The estimator 130 can be implemented in hardware or as a software program that executes in a processor (the same or different processor as the controller 110), which receives the system output 103 at a fixed or variable control cycle sampling interval and uses the new output measurement results and the previous output measurement results to determine the estimated state 121 of the system 120.

[0135] Therefore, by using B&B optimization, the processor reduces the search space by eliminating spaces that do not contain the optimal result. Because the search space is reduced, the processor can accurately determine the optimal solution to control the state of system 120 in a shorter time. Thus, the processor achieves fast processing speed with high precision.

[0136] Figure 1B A block diagram illustrating a predictive controller 110 and a feedback system 120 according to some embodiments is shown. The predictive controller 110 initiates the system 120 such that the estimated state 121 and output 103 of the system 120 follow a command 101. The controller 110 includes a computer, for example, in the form of a single central processing unit (CPU) or multiple CPU processors 151 connected to a memory 152 for storing a system model 102 and constraints 104 on the operation of the system 120. The CPU processor 151 may include a single-core processor, a multi-core processor, a computing cluster, or any number of other configurations. The memory 152 may include random access memory (RAM), read-only memory (ROM), flash memory, or any other suitable memory system.

[0137] Figure 2A A block diagram illustrating a system and method for predictive controller 110, according to some embodiments, for mixed-integer model predictive control (MI-MPC) to compute control signal 111 to achieve a given system current state 121 and control command 101. Specifically, MI-MPC computes control solution 260, such as a solution vector 255 containing a sequence of future optimal discrete and continuous control inputs of the system over the prediction timescale, by solving a constrained mixed-integer optimization problem 250 at each control time step. The objective function, equality, and MIQP data 245 of the discrete and continuous inequality constraints in the optimization problem 250 depend on the dynamic model, system constraints 240, the system current state 121, the control objective, and the control command 101.

[0138] In some implementations, the solution to the inequality-constrained mixed-integer optimization problem 250 uses state and control values ​​at the predicted timescale from a previous control time step 210, which can be read from memory. This concept is referred to as a warm-up or hot-start of the optimization algorithm, and in some implementations, it can reduce the computational effort required by the MI-MPC controller. In a similar manner, the corresponding solution vector 255 can be used to update and store a sequence 260 of optimal state and control values ​​for the next control time step.

[0139] In some implementations, the mixed-integer optimization algorithm is based on a search algorithm, causing the MI-MPC controller to update and store additional mixed-integer programming solution information 260 to reduce the computational workload of the search algorithm in the next control time step. In one implementation, a branch-and-bound optimization method is used to solve the MI-MPC problem at each control time step, and the warm-up information 260 includes data related to nodes in the binary search tree that are part of the solution path from the root node to the leaf node where the optimal integer feasible control solution is found, to improve node selection and variable branching strategies from one control time step to the next.

[0140] Figure 2BA block diagram illustrating a MI-MPC method for solving an optimal control structured mixed-integer optimization problem 250 to compute control signals 111 at each control time step is presented, given the current state 121 and command 101 of system 120. Some implementations are based on a linear dynamic model 263 of the system with linear equality constraints 262, linear continuous inequality constraints 264, linear discrete equality constraints 265, linear terminal inequality constraints 266, and a linear quadratic objective function 261, requiring the solution of a constrained mixed-integer quadratic programming (MIQP) 250 at each control time step. The MIQP data 245 then includes a Hessian matrix and constraint Jacobian matrices 246, along with corresponding gradients and constraint evaluation vectors 247. Generally, the linear discrete equality constraints 265 describe a linear function E of the state and control values. k x k +F k u k The constraint is one of the discrete sets that is equal to the value. In some implementations, the linear discrete equality constraint 265 is a binary equality constraint E. k x k +F k u k For example, the formula includes constraints that a particular control input variable is equal to 0 or 1.

[0141] Figure 3A A schematic diagram illustrating an example of a binary control variable search tree representing a nested tree of search regions for feasible control solutions of integers, according to some implementation methods. Figure 3A A schematic representation of a branch and bound method for implementing a MI-MPC controller is shown in some embodiments, illustrated by a binary search tree 300 in a specific iteration of a mixed-integer optimization algorithm. The main idea of ​​the branch and bound (B&B) method is to sequentially create partitions of the original problem and then attempt to solve these partitions, where each partition corresponds to a specific region of the search space for discrete control variables. In some embodiments, the branch and bound method selects a partition or node and chooses discrete control variables to branch that partition into smaller partitions or search regions, resulting in a nested tree of partitions or search regions.

[0142] For example, partition P1 301 represents a discrete search region that can be divided or branched into two smaller partitions or regions P2 302 and P3 303, i.e., a first region and a second region nested within a common region. The first and second regions are disjoint, i.e., their intersection is empty, P2∩P3=φ307, but together they form the original partition or region P1, i.e., after branching, the union P2∪P3=P1 306 holds. The branch-and-bound method then solves the integer relaxation MPC problem for the first or second partition or region of the search space, obtaining two solutions (local optima) that can be compared with each other and with the currently known upper bound of the optimal objective value. If the performance metrics of the first or second partition or region are not as optimal as the currently known upper bound of the optimal objective value of the MI-MPC problem, then the first or second partition or region can be eliminated. If the first region, the second region, or both regions yield a discrete feasible solution to the MI-MPC problem, then the upper bound can be updated. Then, the branch-and-bound method continues to further partition the region by selecting the remaining regions in the current nested tree of the region.

[0143] While solving each partition may still be challenging, obtaining a local lower bound on the optimal objective value by solving the local relaxations of mixed-integer programming or using duality is quite efficient. If the MI-MPC solver happens to obtain an integer feasible solution while solving the local relaxations, it can use that integer feasible solution to obtain a global upper bound on the mixed-integer controlled solution of the original MI-MPC problem. This helps avoid solving or branching certain partitions that have already been created; that is, these partitions or nodes can be pruned. This general partitioning algorithm idea can be represented as a binary search tree 300, which includes a root node at the top of the tree (e.g., P1 301) and leaf nodes at the bottom of the tree (e.g., P4 304 and P5 305). Furthermore, nodes P2 302 and P3 303 are usually called the direct children of node P1 301, while node P1 301 is called the parent of nodes P2 302 and P3 303. Similarly, nodes P4 (304) and P5 (305) are children of their parent node P2 (302).

[0144] Figure 3B A block diagram is shown for a branch-and-bound mixed integer optimization algorithm that searches for integer feasible optimal control solutions based on nested trees of search regions and corresponding lower / upper bounds, according to some implementations. Figure 3BThe block diagram of the branch-and-bound mixed-integer optimization algorithm shown can be used in some implementations to implement a MI-MPC controller. The branch-and-bound method initializes branch search tree information 310 for mixed-integer quadratic programming (MIQP) at the current control time step based on MIQP data 245, consisting of MIQP matrix 246 and MIQP vector 247. Initialization may additionally utilize branch search tree information and MIQP solution information from previous control time steps 210 to generate a warm-up initialization 310 for the current control time step. The primary objective of the optimization algorithm is to construct lower and upper bounds for the target values ​​of the mixed-integer control solution. In step 311, if the gap between the lower and upper bound values ​​is less than a specific tolerance value, the optimal mixed-integer control solution 255 is found.

[0145] As long as the gap between the lower and upper bounds in step 311 is greater than a specific tolerance value, and the optimization algorithm has not yet reached its maximum execution time, the branch and bound method continues to iteratively search for the mixed-integer optimal control solution 255. Each iteration of the branch and bound method begins by selecting the next node in the tree with possible variable fixing, corresponding to the next region or partition of the integer variable search space, based on a pre-solution branching technique 315. After node selection, the corresponding integer relaxation MPC problem 320 is solved using possible variable fixing, based on a post-solution branching technique.

[0146] If the integer relaxation MPC problem has a feasible solution, then the obtained relaxation control solution provides a lower bound on the target value for that specific region or partition of the integer variable search space. In step 321, if the target is determined to be greater than the currently known upper bound of the target value of the optimal mixed integer control solution, then the selected node 340 is removed from the branch tree. However, in step 321, if the target is determined to be lower than the currently known upper bound, and the relaxation control solution is integer feasible 325, then in step 330 the currently known upper bound and the corresponding mixed integer control solution estimate are updated.

[0147] If the integer relaxation MPC problem has a feasible solution and the objective is below the currently known upper bound 321, but the relaxation control solution is not yet integer feasible, then for the remaining leaf nodes in the branch tree, the global lower bound of the objective can be updated to the minimum objective value 335, and the selected node 340 is removed from the tree. Furthermore, starting from the current node, discrete variables with fractional values ​​are selected for branching according to a specific branching strategy 345, so as to create and append subproblems corresponding to the resulting region or partition of the discrete search space as children of that node in the branch tree 350.

[0148] A crucial step in the branch-and-bound method is creating partitions, specifically selecting which node 315 and which discrete variable to use for branch 345. Some implementations are based on the following: one of the binary control variables with fractional values ​​in the branch-integer relaxation MPC solution. For example, if a particular binary control variable u i,k If ∈{0,1} has a fractional value that is part of the integer relaxed MPC solution, then some implementations add equality constraints u to a subproblem respectively. i,k =0 and add an equality constraint u to another subproblem. i,k =1, to create two partitions for the mixed integer programming. Some implementations are based on a reliability branching strategy 345 for variable selection, which aims to predict future branching behavior based on information from prior branching decisions.

[0149] Some implementations are based on branch-and-bound methods using a depth-first node selection strategy, which can be implemented using a Last-In-First-Out (LIFO) buffer. The next node to be solved is selected as one of the children of the current node, and this process is repeated until a node is pruned—that is, a node that is infeasible, optimal, or dominated by the currently known upper bound—followed by a backtracking process. Conversely, some implementations are based on branch-and-bound methods using a best-first strategy that selects nodes using the current lowest local lower bound. Some implementations combine depth-first and best-first node selection, where a depth-first strategy is used until an integer feasible control solution is found, followed by a best-first strategy in subsequent iterations of the branch-and-bound optimization algorithm. The motivation for the latter is to find an integer feasible control solution as early as possible at the start of the branch-and-bound process (depth-first) to allow for early pruning, followed by a more greedy search for better feasible solutions (best-first).

[0150] The branch and bound method continues to iterate until one or more of the following conditions are met:

[0151] 1. The processor's maximum execution time is reached.

[0152] 2. All nodes in the branch search tree have been removed, so that no new nodes can be selected to solve for convex relaxation or branching.

[0153] 3. The optimal distance between the global lower and upper limits of the objective in a mixed-integer controlled solution is less than the tolerance.

[0154] Figure 4AThe original and dual formulas of convex relaxation, solved as part of a branch-and-bound optimization algorithm in a MI-MPC controller, are illustrated according to some embodiments. When solving a constrained mixed-integer quadratic programming (MIQP), in some embodiments, the original formula 320 of convex relaxation in the branch-and-bound method corresponds to a convex quadratic programming (QP) 400 having a linear quadratic objective function 401, affine inequality constraints 402, and affine equality constraints 403. Some embodiments of this disclosure are based on the understanding that the original optimization variable z = [x, y] between the variables x and y 404 entering the linear quadratic or purely linear objective function 401 can be partitioned. Thus, even if the Hessian matrix H ≥ 0 of the original objective function 401 is positive semi-definite, the matrix... It is also strictly positive definite, and therefore invertible, i.e., Q > 0. Typically, changes in the variable can be used to reform any QP as... Figure 4A The original form.

[0155] Alternatively, in some embodiments, the dual formula of convex relaxation 320 is used in the branch-and-bound method, which corresponds to a dual convex QP 410 with a linear quadratic objective function 411, an affine equality constraint 412, and a Lagrange multiplier positive inequality constraint (μ≥0) 413. The dual optimization variable ν=[μ,λ]415 includes the Lagrange multiplier μ corresponding to the affine inequality constraint 402 and the Lagrange multiplier λ corresponding to the affine equality constraint 403 in the original QP formula 400. Affine function 414 can be used to define a compact notation for the linear quadratic dual objective function ψ(μ,λ)411, which depends on the inverse Q of the positive definite matrix Q>0. -1 .

[0156] Figure 4B The connection 405 between the original formula 400 and the dual formula 410 for each convex relaxation within a branch-and-bound optimization algorithm according to some implementations is illustrated. More specifically, Figure 4BThe relationship between the primal objective function f(z) 401 and the dual objective function ψ(ν) 411 is illustrated. The dual objective function is always concave 430, making the dual QP formula 410 always convex. When solving the convex primal QP 400, the primal objective function is also a convex function 420. The implementation of this disclosure is based on the understanding that if the dual optimization variable ν = [μ,λ] 415 satisfies the affine equality constraint 412 and the positive inequality constraint 413 in the dual QP 410, and if the primal optimization variable z = [x,y] 404 satisfies the affine inequality constraint 402 and the affine equality constraint 403 in the primal QP 400, then the dual objective value of the dual feasible point is less than or equal to the primal objective value of the primal feasible point, i.e., ψ(μ,λ) ≤ f(x,y) 445. Therefore, the dual objective value of the dual feasible point is less than or equal to the dual variable ν. ＊ The dual objective of the optimal value of 431 is 435, and the dual variable is v. ＊ The dual objective 435 of the optimal value of 431 is less than or equal to the original optimization variable z. ＊ The original objective of the optimal value of 421 is 425, that is, ψ(v) ≤ ψ(ν). ＊ )≤f(z ＊ )446.

[0157] When the optimal dual objective value 435 is less than the optimal primal objective value 425, it is generally referred to as weak duality, while when the optimal dual objective value is equal to the optimal primal objective value, it is referred to as strong duality. In the latter case, the objective value of the optimal solution to the dual QP 410 is equal to the objective value of the optimal solution to the primal QP 400. Furthermore, it is well known that strong duality applies to many convex optimization problems if the additional condition, namely the Slater condition, is satisfied.

[0158] Some implementations are based on the understanding that, since ψ(ν) ≤ ψ(ν) in the case of weak duality... ＊ )≤f(z ＊ )446 and in the case of strong duality, ψ(ν)≤ψ(ν) ＊ )=f(z ＊ Therefore, any dual feasible point ν = [μ,λ], that is, the set of values ​​of Lagrange multipliers that satisfy the affine equality constraint 412 and the positive inequality constraint 413 in the dual QP 410, can be used to compute the lower bound of the optimal objective value of the original convex QP 400.

[0159] Some implementations are based on the understanding that when the corresponding node or partition 340 is removed, i.e., when convex relaxation is infeasible or the optimal objective value is greater than the current upper bound (UB) 440, the optimal solution of convex relaxation 320 is not required within the branch and bound optimization algorithm. Additionally, some implementations are based on the understanding that when the optimal objective value is greater than the current upper bound, it is generally easier to find a suboptimal dual solution estimate of ψ(ν)>UB rather than finding the optimal solution of the original formula 400 or dual formula 410 for convex relaxation. Finally, some implementations are based on the understanding that if there is no value z = [x,y] satisfying the affine inequality constraints 402 and affine equality constraints 403 in the original QP 400, then the dual objective function 411 of the dual QP 400 is unbounded when the original QP 400 is infeasible (i.e., the dual objective function 411 is unbounded).

[0160] Figure 4C This illustration demonstrates how, according to some implementations, the suboptimal dual solution estimation, based on its estimation or projection, prematurely terminates the solution process for convex relaxation 320 within a branch-and-bound optimization method. The embodiments of this disclosure are based on the understanding that the suboptimal dual solution estimation ν1 466 can be used to evaluate the dual objective function 411 such that if the corresponding dual objective value 465 satisfies ψ(ν1) > UB, then ψ(ν1) > UB. ＊ If φ(ν1) ≥ ψ(ν1) > UB, then the QP solution process can be terminated early. More specifically, if a suboptimal dual solution estimate is found where the corresponding dual objective value is greater than the current upper limit of 440, this implies that the optimal objective value 435 is also greater than the current upper limit of 440, allowing the corresponding node or partition 340 to be removed within the branch and bound optimization algorithm. Otherwise, the iterative process of the interior point method continues, starting from the suboptimal dual solution estimate, to find the optimal solution for convex relaxation.

[0161] Furthermore, some implementations are based on the understanding that a dual infeasibility estimate ν2 456 might lead to a dual objective value 450 that does not correspond to a point on the concave function 430 defined by the dual objective function 411 of all dual feasible solutions. This means the dual objective value 450 is not necessarily a lower bound of the optimal objective value 435, and therefore the same early termination procedure cannot be applied to such a dual infeasibility estimate. Instead, some implementations perform a computationally efficient projection step 460 for a specific dual infeasibility estimate ν2 456 to obtain a suboptimal dual solution estimate that can be used to evaluate the dual objective function 411. 457, such that if the corresponding dual objective value 455 satisfies Make This allows for the early termination of the QP solution process.

[0162] Figure 5AA block diagram illustrating the process of early termination of the solution to convex relaxation 320 within a branch-and-bound optimization method based on a dual feasible solution estimate, according to some implementations, is provided. This block diagram considers any convex optimization algorithm for solving convex relaxation 320 within the branch-and-bound method, based on an initial solution estimate 510 and by performing one or more convex solver iterations 530. After one or more convex solver iterations 530, a condition can be verified to determine whether a particular solution estimate forms a solution to the convex relaxation 535, in which case a solution to the convex relaxation 540 has been found. Otherwise, if the convex solver termination condition 535 has not been met, another set of convex solver iterations 530 can be performed until a solution 540 is found, indicating that the convex relaxation has no feasible solution, or until the maximum computation time has been reached.

[0163] If the optimal solution 535 for convex relaxation has not yet been found, some implementations evaluate the dual objective function 411 and the dual feasibility constraints 500, including the affine equality constraint 412 and the positive inequality constraint 413 in the dual QP 410. If the dual objective value is greater than the current upper limit 440 and if the solution estimate is dually feasible 505, the solution to the convex relaxation can be terminated early 550, i.e., because the corresponding node or partition is removed 340 in the branch and bound optimization algorithm, so there is no need to compute the optimal solution for the convex relaxation. If the dual objective value is not greater than the current upper limit 440 or the solution estimate is not yet dually feasible 505, another set of convex solver iterations 530 can be executed until a solution is found 540, it is detected that there is no feasible solution for the convex relaxation, the maximum computation time has been reached, or until a dual feasible solution estimate with a dual objective value greater than the current upper limit 505 is computed, causing the convex solver to terminate early 550.

[0164] Figure 5B A block diagram illustrating the early termination process for solving convex relaxations within a branch-and-bound optimization method, based on a projection step that calculates a dual feasible solution estimate, is illustrated according to some embodiments. In the absence of an optimal solution 535 for the convex relaxation, some embodiments evaluate the dual objective function and dual feasibility constraints 500. If the dual objective value is greater than the current upper bound 440 and the solution estimate is dually feasible 505, the solution to the convex relaxation can be terminated early 550, i.e., since the corresponding node or partition is pruned 340, there is no need to calculate the optimal solution for the convex relaxation. If the latter condition 505 is not yet satisfied, the process checks whether the dual objective value ψ is greater than the current upper bound 506. In this case, the projection step 520 can be used to calculate the dual feasible solution estimate; otherwise, another set of convex solver iterations 530 is executed.

[0165] Based on the projection 520 onto the dual feasible solution estimate, some implementations evaluate the dual objective function and dual feasibility constraints 501 of the projection's dual feasible solution estimate. If the new dual objective value ψ +If the solution estimate is greater than the current upper limit and is dually feasible (507), the solution to the convex relaxation can be terminated early (550). Otherwise, another set of convex solver iterations (530) can be executed until a solution is found (540). If it is detected that there is no feasible solution for the convex relaxation, the maximum computation time has been reached, or until a dual feasible solution estimate with a dual objective value greater than the current upper limit is calculated, the convex solver terminates early (550). Therefore, the early termination of the convex solver (550) improves the efficiency of processors configured to determine the optimal solution with low cost (time and storage space) and high accuracy by using branch and bound optimization.

[0166] Some implementations are based on the understanding that evaluating the dual objective function and dual feasibility 500 is generally more computationally cheaper than the projection step 520, which estimates the dual feasible solution, such that the projection step 520 is only performed if the dual objective value ψ of the current solution estimate is greater than the current upper limit 506. Furthermore, some implementations are based on the understanding that if the dual objective value ψ of the solution estimate before the projection step is greater than the upper limit 506, then the dual objective value ψ of the solution estimate after the projection step 520... + It may exceed the current upper limit of 507.

[0167] Figure 6A An example is given of a block sparse interior-point optimization algorithm for solving constrained optimal control structured quadratic programming (OCP-QP) by iteratively solving a relaxed system with necessary optimality conditions, according to some implementations. Furthermore, Figure 6A An example is illustrated of an interior-point optimization algorithm for solving convex relaxations 320 within a branch-and-bound method, which iteratively solves the necessary optimality condition 600 for a smooth system by using a Newton-type solution strategy for the convergence sequence of barrier relaxation parameter values. In some embodiments, within the branch-and-bound optimization algorithm for solving the optimal control structured mixed integer quadratic programming (MIQP) 250, each convex relaxation 320 is a constrained optimal control structured quadratic programming problem (OCP-QP) 610. The OCP-QP optimization problem includes an initial state value condition 612, a system dynamic model leading to linear equality constraints 613, linear inequality constraints 614, and a linear quadratic objective function 611, such that constrained OCP-QP 610 needs to be solved in each convex relaxation 320. The necessary optimality condition 600, sometimes referred to as the Karush-Kuhn-Tucker (KKT) condition, may include one or more stationarity conditions 601, one or more conditions for primal feasibility with respect to equality constraints 602, one or more conditions for primal feasibility with respect to inequality constraints 603, one or more dual feasibility conditions 605, and one or more complementary conditions 604.

[0168] Some implementations are based on the understanding that a relaxed, smooth system is defined by block-structured sparse OCP-QP data 245, containing OCP-QP vectors 247 and OCP-QP matrices 246, in the form of a block-structured sparse Hessian and constraint Jacobian matrix 615, for constraining the necessary optimality condition 600 of the optimal control structured quadratic programming problem (OCP-QP) 610. In some implementations, the interior-point optimization algorithm uses a block-sparse structure utilizing linear algebra operations to reduce the computational complexity of Newtonian-type methods for implementing the smooth system sequence for solving the necessary optimality condition 600.

[0169] exist Figure 6A In the example, the optimization variable z = [x0, u0, x1, u1, ..., x N-1 ,u N-1 ,x N This includes both the state and control input variables of the constrained optimal control structured quadratic programming problem (OCP-QP) 610. The dual variable λ = [λ0, λ1, ..., λ] N [] represents the Lagrange multiplier at each time step k in the control program with respect to equality constraints (e.g., initial state value condition 612 and equality constraint 613 that applies system dynamics). The dual variable μ = [μ0, μ1, ..., μ] N ] represents the Lagrange multiplier at each time step k in the control program with respect to the inequality constraint (e.g., the linear inequality constraint 614), and s = [s0, s1, ..., s2]. N ] represents the relaxation variable corresponding to the inequality constraint. Matrices M and S are diagonal matrices defined by the dual variable μ and the relaxation variable s on the diagonal, respectively, and variable τ represents the parameter value of the barrier-type relaxation of the system under the necessary optimality condition.

[0170] Some implementations use relaxation of the system with a necessary optimality condition different from that using a logarithmic barrier function to obtain one or more smooth nonlinear equations solvable by the Newton-type method in the embodiments of this disclosure. Some implementations are based on the understanding that because the barrier parameter values ​​in the complementary condition 604 are close to zero, i.e., for increasingly smaller barrier parameter values ​​τ>0 and τ→0, the solution to the constrained optimal control structured quadratic programming (OCP-QP) 610 needs to be a solution to the smooth system of the necessary optimality condition 600.

[0171] Figure 6BA block diagram illustrating the initialization steps and iterative process of an interior-point optimization algorithm 600 for solving constrained optimal control problems in the implementation of a predictive control system, according to some embodiments, is provided. The initialization steps may use the optimal or approximate solution from the pre-convex relaxation 620 and / or OCP-QP data 245 to compute initial values ​​of the original optimization variables, Lagrange multipliers (dual variables) with respect to equality and inequality constraints, relaxation variables for inequality constraints, and initial barrier parameter values ​​621.

[0172] Based on the initialization of the optimization variable values ​​621, the iterative process of the interior-point optimization algorithm aims to make the residual value of the first-order necessary optimality condition sufficiently small 626, in which case a convexly relaxed (approximate) optimal solution 320 is found. The iterative process begins by evaluating the residual vector of the first-order necessary optimality condition 625, and then checks whether the norm of the residual vector is sufficiently small relative to the tolerance value 626. A convexly relaxed (approximate) optimal solution 320 is found, and the iterative process terminates if the residual value is sufficiently small; however, if the norm of the residual vector is too large and the number of iterations of the interior-point optimization algorithm has not yet reached its maximum value, the iterative process continues.

[0173] The iterative process of the interior-point optimization algorithm 600 calculates a Newtonian search direction 630 for achieving optimal values ​​of the optimization variables in the constraint convex relaxation of the predictive control system by solving a linearized system (further referred to as a linearized KKT system) of Karush-Kuhn-Tucker (KKT) conditions for a set of first-order optimality conditions. Next, the iterative process calculates a step size 635 in the Newtonian search direction that ensures the positivity of the slack variables and Lagrange multipliers (dual variables) with respect to the inequality constraints. In some implementations, the step size 635 is calculated as the largest positive value smaller than the step size in the Newtonian search direction that ensures the positivity of the slack variables and Lagrange multipliers with respect to the inequality constraints.

[0174] Based on the Newton-type search direction 630 and the calculated step size 635, the iterative process of the interior-point optimization algorithm continues by updating the values ​​of the original optimization variables, the Lagrange multipliers (dual variables), and the relaxation variables 640. Given a new solution estimate for the optimization variables, an update to the barrier parameter values ​​645 can be calculated, allowing the evaluation of a new residual vector for the first-order necessary optimality condition 625. The iterative process of the interior-point optimization algorithm continues until the residual value is sufficiently small 626 and an optimal solution 320 with convex relaxation is found, until the maximum number of interior-point iterations has been reached, or until the interior-point optimization algorithm can be terminated early according to some embodiments of this disclosure.

[0175] Figure 6CA block diagram illustrating a linearized KKT system for predictive control, in which a Newtonian search direction 650 is computed during the iterative process of an interior-point optimization algorithm according to some implementations, to solve the constrained OCP-QP within the branch-and-bound optimization algorithm. The linearized system with KKT conditions can be represented as a block-structured linear equation system 651, which defines the Newtonian search direction Δz of the original optimization variables in the k-th iteration of the interior-point optimization algorithm. k Regarding the Lagrange multiplier (dual variable) Δλ for equality constraints k Regarding the Lagrange multiplier (dual variable) Δμ for inequality constraints k And the Newtonian search direction Δs for slack variables. k The evaluation of the stationarity condition in the k-th iteration of the interior point optimization algorithm for the right-hand vector of the linear KKT system 651 is as follows: 652; 653; 654; 655; 655; 656; 657; 658; 659; 650; 651; 652; 653; 654; 655; 655; 655.

[0176] The implementation of this disclosure is based on the following understanding: for continuous iteration of the interior point optimization algorithm, for non-zero obstacle parameter values ​​τ k >0, the relaxed complementary condition 655 corresponds more closely to the increasingly smaller values ​​of the barrier parameter τ. k →0 is an exact complementary condition.

[0177] Figure 6D A block diagram illustrating a sequence of two subsystems in an interior-point optimization algorithm for solving constrained OCP-QP, according to some implementations, is shown. Some implementations are based on the understanding that a linearized KKT system can be decomposed into a symmetric linear subsystem 661 to compute the Newtonian search direction Δz of the original optimization variables. k Regarding the Lagrange multiplier (dual variable) Δλ for equality constraints k And the Newton-type search direction Δμ with respect to the Lagrange multipliers (dual variables) of the inequality constraints. k In the k-th iteration of the interior-point optimization algorithm, after solving the latter symmetric linear subsystem 661, the diagonal linear system 662 can be solved to calculate the Newtonian search direction Δs of the slack variables. k .

[0178] The right-hand vector of the first symmetric linear KKT system 661 consists of the following: an evaluation of the stationarity condition 652, an evaluation of the original feasibility of the equality constraints 653, and the Newtonian search direction Δs due to the slack variables. kThis leads to a modified variant 663 of the original feasibility assessment 654 relative to the inequality constraints. Similarly, the right-hand vector of the diagonal linear system 662 contains, in addition to the Newton-type search direction Δμ due to the Lagrange multipliers (dual variables) of the inequality constraints. k Evaluation of the relaxation complementarity condition 655 caused by the contribution.

[0179] Some implementations are based on the understanding that the first symmetric linear KKT system 661 can be solved based on iterative linear algebra routines or by direct factorization or decomposition of the symmetric, sparse, and block-structured KKT matrices on the left side of the system. Furthermore, some implementations are based on the understanding that for the k-th iteration of the interior-point optimization algorithm, due to the Lagrange multipliers (dual variables) μ on the diagonal with respect to the inequality constraint... k Defined matrix M k The diagonal sparse structure allows each equation of the diagonal linear system 662 to be solved independently based on a single scalar division.

[0180] Figure 7A A block diagram illustrating the process of calculating the projection step of the suboptimal dual solution estimate 520 within the branch-and-bound optimization algorithm, according to some implementations, can be used to prematurely terminate the convex solver. This is for the Lagrange multipliers (μ) that are nearly dualally feasible or dualally infeasible. k ,λ k Given the current set of values ​​701, compute the suboptimal dual solution estimate 520. First, compute the update (Δμ) of the Lagrange multiplier values. k ,Δλ k )710. Then, the process verifies for the updated value (μ) k +Δμ k ,λ k +Δλ k Whether the dual feasibility constraint 705 is satisfied, that is, whether the solution estimate satisfies the affine equality constraint 412 and the positive inequality constraint 413 in the dual QP formula 410 of the convex relaxation 320. If the solution estimate updated after the projection step is still not dual feasible 705, the projection step is unsuccessful 715 and the process continues by performing one or more convex solver iterations 530 until an optimal solution 540 is found, it is detected that there is no feasible solution for the convex relaxation, or until the maximum computation time has been reached.

[0181] If the value after projection is dually feasible (705), then the current value of the Lagrange multiplier (720) is updated, and the process checks whether these updated values ​​can be used to prematurely terminate the convex solver (725) within the branch and bound optimization algorithm. The implementation of this disclosure is based on the understanding that if the solution estimate updated after the projection step is dually feasible and the dual objective is greater than the current upper bound (507), then the updated value of the Lagrange multiplier (720) can be used to prematurely terminate the convex solver (550).

[0182] Figure 7B A block diagram illustrating the process of calculating the projection step of the suboptimal dual solution estimate 520 by solving an equality-constrained linear quadratic optimization problem to enforce affine equality constraints 412 in the dual QP formula 410 of convex relaxation 20 is illustrated according to some embodiments. Some embodiments are based on the understanding that the projection step (Δμ) of the Lagrange multipliers can be calculated by solving an equality-constrained quadratic programming (QP) optimization problem 730. k ,Δλ k )710, where the quadratic objective makes the updated size (Δμ) k ,Δλ k )731 is minimized, and equality constraint 732 enforces affine constraint 412 in dual QP 410. To this end, the projection step computes an update to the suboptimal dual solution estimate 520 by solving an equality-constrained optimization problem that enforces the real dual feasibility constraint and minimizes the update norm of the Lagrange multiplier variables.

[0183] Some implementations are based on the understanding that the equality constraint QP of projection step 730 can be effectively solved by using structured leverage direct or iterative linear algebra routines with computational costs significantly lower than solving convex relaxation 320 to solve block-structured linear systems. However, the projection step (Δμ) based on equality constraint QP... k ,Δλ k )730 does not directly enforce positive inequality constraints 413 in the dual QP 410, making the updated values ​​of the Lagrange multipliers potentially not always dually feasible 705, and the projection step potentially not always successful 715. In some implementations, the projection step 730 computes an update to the estimate of the suboptimal dual solution by solving a block-structured linear system whose solutions form the optimal solution to the equality-constrained optimization problem.

[0184] In some embodiments of this disclosure, the projection step 710 is computed by performing one or more iterations of the process of solving the equality constraint QP 730 until a dual feasible solution estimate 520 that can be used to prematurely terminate the convex solver 725 is computed, or until the maximum computation time has been reached, allowing the process to continue by performing one or more convex solver iterations 530. Embodiments of this disclosure are based on the understanding that matrix factorization can be reused in subsequent solutions to the equality constraint QP 730 to reduce the computational cost of computed projection steps for prematurely terminating the dual feasible solution estimate of the convex solver.

[0185] Figure 7C A block diagram illustrates the process of calculating the projection step of the dual feasible solution estimate 520 by solving an equality-constrained linear quadratic optimization problem using alternative scaling of the optimization variables in the objective of an equality-constrained QP that effectively enforces the updated Lagrange multiplier values, according to some embodiments. Furthermore, the equality-constrained linear quadratic optimization problem is solved by attempting to enforce real additive positive inequality constraints 413 in the dual QP formula 410 of convex relaxation 320, while simultaneously using alternative scaling of the optimization variables in the objective that effectively enforces real additive affine equality constraints 412. Some embodiments are based on the understanding that the projection step (Δμ) of the Lagrange multipliers can be calculated by solving the equality-constrained quadratic programming (QP) optimization problem 740. k ,Δλ k )710, including the update Δμ of the dual variable as the optimization variable in equality constraint QP 740. k ,Δλ k and the update of the original variable Δx k .

[0186] Figure 7C The projection step also enforces one or more additional optimality conditions; and minimizes the weighted norm of the updates to the original optimization variables and the Lagrange multiplier variables based on the positive definite weighting matrix in the objective function. The positive definite weighting matrix in the objective of the projection step corresponds to the Hessian matrix used for the weighting of the updates to the original optimization variables, the diagonal matrix used for the small positive weighting of the updates to the Lagrange multipliers corresponding to the equality constraints in the convex relaxation, and the positive weighting corresponding to the ratio between the relaxation variables and the Lagrange multiplier values ​​corresponding to each inequality constraint after one or more iterations in the convex solution process.

[0187] The quadratic objective makes the updated size (Δx) k ,Δμ k ,Δλ k Minimize 741 using the positive definite Hessian matrix from the original QP formula 400 derived from convex relaxation. Right now Use the target item Small positive values ​​∈ dual >0 minimizes the update step size of the Lagrange multipliers corresponding to the equality constraints. The objective term is obtained using... Positive definite weighted matrix This minimizes the update step size of the Lagrange multipliers corresponding to the inequality constraints. Finally, equality constraint 742 enforces affine constraint 412 in the dual QP 410, as well as additional constraints designed to enforce that the updated solution estimate still closely satisfies the necessary optimality condition corresponding to convex relaxation (i.e., the Karush-Kuhn-Tucker (KKT) condition).

[0188] Some implementations use an interior-point optimization algorithm 600 to solve each convex relaxation within the branch-and-bound optimization algorithm, such that in the k-th iteration of the interior-point optimization method (IPM) 600 for solving constrained convex relaxation problems, for the values ​​of the relaxation variables... With the value of Lagrange multipliers The ratio between For each corresponding inequality constraint, the positive definite weighting matrix in the quadratic objective 741 of the equality constraint QP in projection step 740. It can be calculated as a diagonal matrix with diagonal elements.

[0189] Some implementations are based on the understanding that if an inequality constraint makes equality true for the optimal value of the optimization variable, then the inequality constraint is considered "active" at the optimal solution of a convex relaxation. If an inequality constraint makes inequality true for the optimal value of the optimization variable, then the inequality constraint is considered "inactive" at the optimal solution of a convex relaxation.

[0190] Some implementations are based on the understanding that in a series of subsequent iterations of the interior-point optimization algorithm, for the slack variable s i >0 is non-zero and the corresponding Lagrange multiplier value remains positive but gets smaller and smaller and approaches zero (i.e., μ). i →0) However, for inequality constraints that are inactive at the optimal solution, the value of w grows larger and larger until it reaches infinity, i.e., w i →∞. Even if the current value of the Lagrange multiplier is... 701 can be relatively close to zero. For inequality constraints that are inactive at the optimal solution, the objective term... Weighted matrix in 741 diagonal The value of increases progressively with subsequent iterations of the interior point solution process, causing the corresponding update from projection step 740 to... It may be very small and may satisfy the positive constraint. 705. If projection 740 successfully makes the updated value of the Lagrange multiplier positive and satisfies the duality feasibility constraint 705, then the solution estimate can be used to check for early termination 725 of the convex solver to reduce the computational cost of the branch-and-bound optimization algorithm in the predictive controller.

[0191] Some implementations are based on the understanding that in a series of subsequent iterations of the interior-point optimization algorithm, for the Lagrange multiplier value μ... i >0 is non-zero and the corresponding value of the slack variable remains positive but becomes smaller and smaller to approach zero (i.e., s i →0) However, for the inequality constraints active at the optimal solution, the value of w remains positive, but it becomes smaller and smaller, tending towards zero, i.e., w i →0. For inequality constraints active at the optimal solution, the objective term... Weighted matrix in 741 diagonal The value of becomes increasingly smaller with each subsequent iteration of the interior point solution process, and allows for corresponding updates from projection step 740. Relatively large because of the Lagrange multipliers The current value 701 is positive and relatively large, making it possible to satisfy the positive constraint. 705, which allows the convex solver to terminate earlier by 725.

[0192] Figure 7D Pseudocode is provided for a procedure 750, according to some implementations, for early termination of convex relaxation solutions within a branch-and-bound optimization algorithm in a predictive controller, based on evaluation of the dual objective and projection steps on estimation of dual feasible solutions. The solution process for each convex relaxation can begin from the original optimization variable z. 0 Lagrange multipliers μ 0 and λ 0 and slack variables s 0The process begins with a fixed set of warm-start solution estimates or values, and a current upper bound value 751 is given during the branch-and-bound optimization process. The solution process continues as long as the residual value is not sufficiently small (626), in which case the process first evaluates the dual objective function and dual feasibility constraints (500). If the dual objective value is greater than the current upper bound (440), and if the solution estimate is dually feasible (505), the solution to the convex relaxation can be terminated early (550), i.e., the optimal solution to the convex relaxation does not need to be computed because the corresponding nodes or partitions are pruned (340). If the latter case (505) is not yet satisfied, the process checks whether the dual objective value ψ is greater than the current upper bound (506), in which case the projection step can be used to compute the dual feasible solution estimate (710); otherwise, another set of convex solver iterations (530) is executed. In some implementations, the latter convex solver iterations (530) correspond to iterations of the interior point optimization algorithm (600) that solves each constraint convex relaxation within the branch-and-bound method.

[0193] Based on the calculated projection step 710, updated values ​​can be calculated for the Lagrange multipliers 752 that evaluate the dual objective function and dual feasibility constraint 501. If the updated value after the projection step is dually feasible 705, then the current value 720 of the Lagrange multipliers is updated. Furthermore, if the new dual objective value ψ... + If the solution estimate is greater than the current upper limit and is dually feasible (507), the solution to the convex relaxation can be terminated early (550). Otherwise, another set of convex solver iterations can be executed (530) until a solution is found (540). If no feasible solution is found for the convex relaxation, the maximum computation time has been reached, or the convex solver terminates early (550) until a dual feasible solution estimate with a dual objective value greater than the current upper limit is calculated.

[0194] Figure 8 This illustrates the equivalence of the equality-constrained QP solution and the solution of the simplified linear system in the computational projection step according to some implementations, followed by the independent computation of the equations for the approximately linearized KKT system with dual variable updates. Furthermore, Figure 8 The equations of the approximately linearized KKT system 810 are shown. Solving the equations of the approximately linearized KKT system 810 to compute the projection step 710 as the optimal solution to the equality-constrained quadratic programming (QP) optimization problem 740 includes updating Δμ, which is the dual variable of the optimization variable in the equality-constrained QP 740. k ,Δλ k and the update of the original variable Δx k Some implementations are based on the understanding that the optimal solution to the equality constraint QP 740 is equivalent to the solution to the linear equation system 810 of 805.

[0195] Some implementations are based on the understanding that the linear system 810, which needs to be solved to compute the equations for the projection step 740 in order to compute the dual feasible solution estimate for early termination, has a specific block sparse structure whose matrix is ​​very close to the symmetric matrix in the linearized KKT system 661 in the interior-point optimization algorithm 600. However, in the iteration of the interior-point optimization algorithm 600, the right-hand side 811 of the linear system 810 computed in the projection step 710 is different from the right-hand side of the linearized KKT system 661 computed in the Newton-type search direction. The right-hand side 811 of the projection step is zero, except for the stationarity condition 652 concerning the optimization variable y 404 that does not enter the original objective function 401 or enters the original objective function 401 linearly, which corresponds to the dual feasibility equality constraint 412 in the dual QP 410. Some implementations are based on the understanding that duality is enforced on the right side of the projection step 811, while aiming to preserve the ability to continue progressing until the current iteration of the interior point optimization algorithm satisfies the necessary optimality condition corresponding to convex relaxation, namely, the Karush-Kuhn-Tucker (KKT) condition.

[0196] The Hessian matrix H≥0 of the original objective function 401 is positive semi-definite, ∈ dual >0 is a small positive penalty weight value used for the Lagrange multipliers corresponding to equality constraints, while This is the positive definite penalty weighting matrix used for the Lagrange multipliers corresponding to the inequality constraints. Matrices F and G correspond to the constraint Jacobian matrices for the equality and inequality constraints, respectively. In some implementations, in the k-th iteration of the interior-point optimization method (IPM) 600 for solving constraint convex relaxation in the branch-and-bound optimization algorithm, for the values ​​of the relaxation variables... With the value of Lagrange multipliers The ratio between For each corresponding inequality constraint, the weighting matrix It is computed as a diagonal matrix with diagonal elements.

[0197] Some implementations are based on the understanding that the approximate linearized KKT system in projection step 810 is equivalent to effectively solving the simplified linearized KKT system in step 815 to calculate the Newtonian search direction Δz of the original optimization variables 821, followed by the projection step 822 of calculating the Lagrange multipliers. The block-structured linear system composed of linearized systems with KKT conditions is solved by solving the simplified linear system to calculate the update of the original optimization variables, and then solving the diagonal linear system to calculate the update of the Lagrange multiplier variables. Some implementations are based on the understanding that the matrix in the simplified linearized KKT system 821...

[0198]

[0199] It is positive definite, allowing the Cholesky factorization to be computed over the matrix, thus solving the simplified linearized KKT system 821 in a computationally efficient manner. The projection step 822 of the Lagrange multipliers corresponding to the equality constraints can be independently computed as follows: In the k-th iteration of the interior point optimization method (IPM) 600 for constrained convex relaxation in the branch and bound optimization algorithm, in the weighting matrix... When computed as a diagonal matrix, the projection step 822 of the Lagrange multipliers corresponding to the inequality constraints can be independently computed as follows:

[0200] Figure 9A An example is illustrated based on a positive definite approximation KKT matrix M according to some implementations. k Forward Cholesky factorization is used to solve the simplified approximate linearized KKT system 821 to calculate the Newton-type search direction Δz of the original variables. k In order to calculate the projection step Δμ of the Lagrange multipliers k ,Δλ k A flowchart of the process. The initial step is potentially given by the positive definite approximation KKT matrix M in the earlier iterations of the linear subsystem 821 of the IPM. k In the current iteration of the interior-point optimization method (IPM), the positive definite approximation KKT matrix M in the simplified linear subsystem 821 is computed or updated. k 900. Given Cholesky factorization information from previous IPM iterations, the forward Cholesky factorization or decomposition can be computed or updated. (where matrix L is a lower triangular matrix) 905, and then use it to solve the simplified approximate linearized KKT system 910.

[0201] In some embodiments of this disclosure, as part of the projection step 820, which can be used to terminate the convex solver early, the lower triangular system is first calculated via forward substitution. The solution is 911, and then the upper triangular system is calculated by backward substitution. The solution 912 yields the solution vector of the simplified linearized KKT system. 913, based on forward Cholesky factorization or decomposition 905 performs a simplified approximate linearized KKT system The solution is 910. Therefore, given the Cholesky factorization from a prior projection step or a prior iteration from a convex solution process, the Cholesky factorization is computed or updated to compute the solution of a simplified linear system with a positive definite KKT matrix.

[0202] Some implementations are based on the understanding that the matrix of the simplified approximate linearized KKT subsystem 821, which is defined to compute the Newton-type search direction, is positive definite and exhibits a block tridiagonal sparse structure due to the block-structured sparsity of the constrained optimal control problem, making it possible to solve the simplified approximate linearized KKT system using block tridiagonal Cholesky factorization. In some implementations, the KKT matrix of the simplified linear system has a block tridiagonal sparse structure due to the block-structured sparsity of the mixed-integer optimal control optimization problem, and the block tridiagonal Cholesky factorization is computed or updated to compute the solution of the simplified linear system with a positive definite block tridiagonal KKT matrix. In other implementations of this disclosure, dense, banded, or sparse Cholesky factorization can be used to reduce the computational cost of solving the structured linear system in the projection step of the early termination process.

[0203] Some implementations are based on the understanding that the positive definite weighting matrix corresponding to the inequality constraints in the transition from one IPM iteration to the next IPM iteration... In the case of low-rank updates, the low-rank updates of the KKT matrix and its Cholesky factorization can be computed at a considerably low computational cost, compared to a completely new evaluation of the KKT matrix and its Cholesky factorization.

[0204] Figure 9B The diagram illustrates a method based on a positive definite approximation KKT matrix M according to some implementations. k The inverse Cholesky factorization is used to solve the simplified approximate linearized KKT system 821 to calculate the Newton-type search direction Δz of the original variables. k A flowchart of the process. Potentially given information from prior IPM iterations regarding the reverse Cholesky factorization, the reverse Cholesky factorization or decomposition can be computed or updated. (where matrix R is an upper triangular matrix) 915, and then use it to solve the simplified approximate linearized KKT system 920.

[0205] In some implementations, as part of the projection step 820, which can be used to terminate the convex solver early, the upper triangular system is first calculated via backward substitution. The solution is 921, and then the lower triangular system is calculated through forward substitution. The solution 922 yields the solution vector of the simplified linearized KKT system. 923, based on inverse Cholesky factorization or decomposition 915 Perform the solution of the simplified approximate linearized KKT system. 920.

[0206] Figure 9C This illustration shows a block-structured sparse pattern 930 for a linearized system of Newtonian search directions for calculating the original and dual optimization variables in projection step 710, according to some embodiments. Due to the separable, stage-by-stage structure of the objective in the constrained optimal control structure optimization problem, the objective's Hessian matrix H can be diagonal or block-diagonal 926. Due to the stage-by-stage coupling of state variables at subsequent time intervals in the constrained OCP, the equality constraint Jacobian matrix F 931 and its transpose... 927 typically exhibits a block-diagonal sparse structure. Due to the separable, stage-by-stage structure of inequality constraints on state and / or control input variables in constrained OCP, the inequality constraint Jacobian matrix G 933 and its transpose... 928 can be diagonal or block diagonal.

[0207] The block matrices on the diagonal of the block-structured KKT matrix 925 (such as the block matrix corresponding to equality constraint 932 and the block matrix corresponding to inequality constraint 934) are diagonal matrices. More specifically, the diagonal elements corresponding to equality constraint 932 are equal to –∈ dual , where ∈ dual Values ​​greater than 0 are relatively small. The diagonal elements of 934, corresponding to the inequality constraints, depend on each specific inequality constraint and can be equal to their corresponding... The values ​​can be defined between the lower and upper limits.

[0208] Figure 9D An example is given of the Newton-type search direction Δz for the original variable in the projection calculation step 710 according to some embodiments. kA schematic diagram of the block-structured sparse pattern of the positive (semi-)definite KKT matrix 935 in the simplified approximate linearized KKT system 821 is shown in Figure 940. Some embodiments are based on the understanding that the block tridiagonal sparse structure 941 of the KKT matrix can be used to solve simplified linear systems based on block tridiagonal matrix factorization (e.g., block tridiagonal forward or backward Cholesky factorization) to reduce the computational cost of the projection step in the early termination process. Other embodiments of this disclosure are based on the understanding that the striped sparse structure 941 of the KKT matrix can be used to solve simplified linear systems based on striped and / or ordinary sparse matrix factorization to reduce the computational cost of the projection step in the early termination process in predictive control systems.

[0209] Figure 9E An example is illustrated based on a block tridiagonal positive definite approximation KKT matrix M according to some implementations. k 950 block tridiagonal forward Cholesky factorization 955 solution simplifies approximate linearization of KKT system 821 to calculate Newton-type search direction Δz of original variables k In order to calculate the projection step Δμ of the Lagrange multipliers. k ,Δλ k A flowchart of the process. Given the block tridiagonal forward Cholesky factorization information from previous IPM iterations or from the projection steps in previous iterations of the convex relaxation solution process, the block tridiagonal block forward Cholesky factorization or decomposition can be computed or updated. (where matrix L is the lower block diagonal matrix) 955, and then use it to solve the simplified block tridiagonal linearized KKT system 960.

[0210] In some implementations, the next-diagonal system is calculated first via block forward substitution. The solution is 961, and then the upper block two-diagonal system is calculated by block backward substitution. The solution 962 yields the solution vector of the simplified block tridiagonal linearized KKT system. 963, based on block tridiagonal forward Cholesky factorization or decomposition 955 performs simplified block tridiagonal linearization of KKT system Solve for 960.

[0211] Figure 9F An example is illustrated based on a block tridiagonal positive definite approximation KKT matrix M according to some implementations. k 950 block tridiagonal inverse Cholesky factorization 965, solving simplified approximate linearized KKT system 821 to calculate Newton-type search direction Δz of the original variables. kA flowchart of the process. Given block tridiagonal inverse Cholesky factorization information from previous IPM iterations or from projection steps in previous iterations of the convex relaxation solution process, the block tridiagonal inverse Cholesky factorization or decomposition can be computed or updated. (where matrix R is the upper block diagonal matrix)965, and then use it to solve the simplified block tridiagonal linearized KKT system970.

[0212] In some implementations, the upper block bidiagonal system is calculated first via block backward substitution. The solution is 971, and then the next block two-diagonal system is calculated by block forward substitution. The solution 972 yields the solution vector of the simplified block tridiagonal linearized KKT system. 973 is based on block tridiagonal inverse Cholesky factorization or decomposition 965 is used to perform simplified block tridiagonal linearization of the KKT system. Solve for 970.

[0213] Some implementations are based on the understanding that it is only necessary to compute or update one or more blocks of the inverse block tridiagonal Cholesky factorization 965, which correspond to blocks containing control processes near the start of the constrained optimal control optimization problem in the predictive control system. The inequality constraints on the updated values ​​of the blocks result in a significant reduction in computational cost. Other implementations are based on the understanding that only one or more blocks of the forward block tridiagonal Cholesky factorization 955 need to be computed or updated, the one or more blocks corresponding to the end of the control process in the constrained optimal control optimization problem of the predictive control system. The block of inequality constraints for updating values ​​leads to a significant reduction in computational cost.

[0214] Figure 10A A block diagram illustrates a warm-start initialization process for calculating solution estimates based on approximate or optimal solutions of convex relaxations in parent nodes, according to some implementations. Furthermore, the warm-start initialization process calculates solution estimates for convex relaxations within the branch-and-bound optimization algorithm. This process first reads the solution information 1000 of the convex relaxations in the parent problem of the branch-and-bound search tree from memory, i.e., the optimal or approximate values ​​(z) of the original variables, dual variables, and relaxation variables. 0 ,λ 0 ,μ 0 ,s 0Some implementations are based on the understanding that the dual feasible solution estimate of the parent problem additionally forms the dual feasible solution estimates of one or more of its children in the branch-and-bound search tree, such that this dual feasible solution estimate can be used for early termination to reduce the computational cost of solving convex relaxations in one or more child nodes. In another implementation, the projection step computes a suboptimal dual solution estimate 520 based on the local optimum to perform an early termination process for one or more other regions within the nested region tree of the B&B optimization.

[0215] In some implementations, the warm-start initialization process uses a positive scaling factor κ >> 1 to scale the initial values ​​of the slack variables upwards. Initial values ​​of and / or dual variables 1010, and calculate obstacle parameters. The resulting new value is 1015. Some implementations use an additional saturation function to ensure the lower and upper limits of the initial barrier parameter value, i.e., τ. min ≤τ 0 ≤τ max , where τ min and τ max These represent the expected minimum and maximum values ​​of the barrier parameters, respectively. A set of values ​​obtained can be used. As an initial estimate for the interior point optimization algorithm, it is used to solve the convex relaxation 1020 at the current node in the branch and bound method.

[0216] Some implementations are based on the understanding that upscaling 1010 of one or more values ​​of the slack variable leads to a constraint deactivation process designed to avoid any slow convergence of the interior-point optimization algorithm due to errors in treating inactive constraints as active constraints during warm-up initialization. Some implementations are based on the understanding that upscaling 1010 of one or more values ​​of the dual variable (Lagrange multiplier) leads to a constraint activation process designed to avoid any slow convergence of the interior-point optimization algorithm due to errors in treating active constraints as inactive constraints during warm-up initialization. More specifically, when the scaling parameter value κ >> 1 is chosen to be sufficiently large but not excessively large, some or all of the inequality constraints can be moved away from the non-smoothing of the complementarity condition, which states that for each inequality constraint in the optimal solution of the convex relaxation, the slack variable or Lagrange multiplier value should be equal to zero. Furthermore, the initial estimate 1015 of the barrier parameter increases the set of relaxed and thus smoothed complementarity conditions.

[0217] Figure 10BA block diagram illustrating a smooth relaxation-based warm-start initialization process for solving one or more interior-point optimization algorithms for convex relaxation in a predictive control system, according to some embodiments, is provided. Some embodiments are based on the understanding that when the norm of the residual vector is sufficiently small and / or when the barrier parameter value is sufficiently small, an approximate solution of the barrier-type relaxation, constraining the optimality conditions of the convex relaxation, can be stored in memory for the parent node in the branch-bound search tree. The first step of the smooth relaxation-based warm-start initialization for the convex solution process is then to read the approximate relaxation solution information (z) of the convex relaxation in the parent problem of the branch-bound search tree from memory. 0 ,λ 0 ,μ 0 ,s 0 ,τ 0 )1030. Next, we obtained a set of values ​​(z) 0 ,λ 0 ,μ 0 ,s 0 ,τ 0 It can be used as an initial estimate for interior point optimization algorithms to solve for the convex relaxation 1020 at the current node in the branch and bound method.

[0218] Some implementations are based on the understanding that the performance of a warm-start initialization process based on smooth relaxation is highly dependent on the decision rules to store approximate solutions of the constraint convex relaxation optimality conditions of the barrier relaxation in the parent nodes of the branch-and-bound search tree. For example, if the norm of the residual vector and the barrier parameter values ​​are both below a certain threshold, approximate and relaxed solution information can be stored in a specific iteration of the interior-point optimization algorithm.

[0219] Figure 11 A block diagram illustrating an early termination process for solving convex relaxations within a branch-and-bound algorithm, according to some implementations, is provided when the optimal objective value is greater than the current upper limit or when no feasible solution exists for the convex relaxation. Furthermore, Figure 11The projection step based on the estimation or computation of the dual feasible solution 1100 is illustrated, along with an early termination process 750 that reduces the computational cost of solving convex relaxation within the branch-and-bound optimization algorithm in the predictive controller, based on the evaluation of the dual objective value 1110. The convex relaxation solution process 550 can be terminated early if the dual objective value is greater than the current upper bound and if the solution estimate is dually feasible 1115; that is, the optimal solution to the convex relaxation does not need to be computed because the corresponding node or partition is pruned 340. Some implementations are based on the following: satisfying the latter condition 1120 if the optimal objective value is greater than the current upper bound (UB) in the branch-and-bound optimization algorithm 1130, or if no feasible solution to the constrained convex relaxation problem exists 1140. More specifically, this implementation is based on the following realization: when the original QP 400 is not feasible, the dual objective function 411 of the dual QP 410 is unbounded, that is, if there is no value z = [x,y] in the original QP 400 that satisfies the affine inequality constraint 402 and the affine equality constraint 403, then the dual objective function 411 is unbounded and therefore greater than any upper bound value in the branch and bound method.

[0220] Figure 12A A schematic diagram of a vehicle 1201 is illustrated, incorporating a predictive controller 1202 employing principles from several implementations. As used herein, vehicle 1201 can be any type of wheeled vehicle, such as a passenger car, bus, or taxi. Furthermore, vehicle 1201 can be an autonomous or semi-autonomous vehicle. For example, some implementations control the movement of vehicle 1201. Examples of movement include lateral movement of the vehicle controlled by a steering system 1203 of vehicle 1201. In one implementation, the steering system 1203 is controlled by controller 1202. Additionally or alternatively, the steering system 1203 may be controlled by the driver of vehicle 1201.

[0221] The vehicle may also include an engine 1206 that can be controlled by the controller 1202 or other components of the vehicle 1201. The vehicle may also include one or more sensors 1204 for sensing the surrounding environment. Examples of sensors 1204 include rangefinders, radar, lidar, and cameras. The vehicle 1201 may also include one or more sensors 1205 for sensing its current motion and internal state. Examples of sensors 1205 include Global Positioning System (GPS), accelerometers, inertial measurement units, gyroscopes, shaft rotation sensors, torque sensors, deflection sensors, pressure sensors, and flow sensors. The sensors provide information to the controller 1202. The vehicle may be equipped with a transceiver 1206, enabling the controller 1202 to communicate via wired or wireless communication channels.

[0222] Figure 12BA schematic diagram illustrating the interaction between a controller 1202 (i.e., a predictive controller) and other controllers 1220 of a vehicle 1201 according to some embodiments is shown. For example, in some embodiments, the controller 1220 of the vehicle 1201 is a steering 1225 and a brake / throttle controller 1230 that respectively control the rotation and acceleration of the vehicle 1201. In this case, the predictive controller 1202 outputs control inputs to controllers 1225 and 1230 to control the state of the vehicle 1201. The controller 1220 may also include a higher-level controller, such as a lane-keeping assist controller 1235, that further processes the control inputs of the predictive controller 1202. In both cases, the controller 1220 uses the output of the predictive controller 1202 to control at least one actuator of the vehicle 1201 (such as the steering wheel and / or brakes of the vehicle 1201) to control the movement of the vehicle 1201. Furthermore, the predictive controller 1202 determines the inputs of the vehicle 1201 based on the mixed-integer control solution, wherein the inputs of the vehicle 1201 include the acceleration of the vehicle 1201, the engine torque, braking torque and steering angle of the vehicle 1201, and one or a combination of discrete optimization variables, in order to model one or a combination of discrete control decisions, switching, gear shifting and obstacle avoidance constraints in system dynamics.

[0223] Figure 12C Schematic diagrams illustrating path and / or motion planning methods for controlled vehicles, employing principles from some implementation methods, are provided. Furthermore, Figure 12C A schematic diagram illustrates an autonomous or semi-autonomous controlled vehicle 1250 for which a dynamically feasible and generally optimal trajectory 1255 can be calculated using embodiments of the present disclosure. The generated trajectory is designed to keep the vehicle within a specific road boundary 1252 and to avoid other controlled and / or uncontrolled vehicles, i.e., those vehicles that are obstacles 1251 to the particular controlled vehicle 1250. In some embodiments, each obstacle 1251 may be represented by one or more inequality constraints in the time or space formula of a mixed-integer programming problem, including one or more additional discrete variables for each obstacle. For example, based on an embodiment configured to implement a mixed-integer model predictive controller, the autonomous or semi-autonomous controlled vehicle 1250 may make discrete decisions in real time, such as overtaking on the left or right, or instead, staying behind another vehicle in the current lane of road 1252, while simultaneously making continuous decisions in real time, such as speed, acceleration, or steering inputs, to control the movement of the vehicle 1250.

[0224] Figure 12D Exemplary traffic scenarios are illustrated based on single-vehicle or multi-vehicle decision-making modules according to some implementation methods. Figure 12DThe following scenario is depicted: There is one or more controlled vehicles (referred to as ego vehicle 1271), and traffic consisting of other vehicles similarly shown in 1272, such as lane 1273 marked L6, stop line 1274 marked S1, and intersection 1275 marked I3. For a vehicle at location 1261 with a final destination 1262, the routing module provides a road sequence indicated by arrow 1263 and a turning sequence indicated by arrow 1264. However, it should be noted that the road sequence 1263 and the turning sequence 1264 do not themselves specify the vehicle's trajectory or path. There are many discrete decisions to be made, such as which lane the vehicle should travel in, whether the vehicle should change lanes or remain in the current lane, whether the vehicle should begin to decelerate to stop at the stop line, whether the vehicle is allowed to cross the intersection, and so on. Furthermore, there are many continuous decisions to be made, such as the time series of positions and directions the vehicle should reach during its journey from its initial point to its destination. These decisions depend heavily on the current traffic conditions at the moment the vehicle arrives at the corresponding location. Due to the uncertainty of traffic movement and the uncertainty of when the vehicle will arrive at that location, the routing module is typically unaware of the current traffic conditions. In some embodiments of this disclosure, a motion plan can be computed for one or more controlled self-vehicles 1271 that may have communication to allow coordination between vehicles (V2V) and / or between intelligent infrastructure systems and vehicles (V2X).

[0225] Figure 13A and Figure 13B This is a schematic diagram illustrating the formulation of a spacecraft mixed-integer predictive control problem based on the principles of some embodiments of this disclosure. More specifically, Figure 13A and Figure 13B An example is illustrated of a spacecraft 1302 equipped with multiple actuators such as thrusters 1350 and momentum exchange devices 1351. Examples of momentum exchange device types include reaction wheels (RWs) and gyroscopes. Spacecraft 1302 is a vehicle, container, or machine designed to fly in space, whose operation changes quantities such as position, velocity, attitude, or orientation of spacecraft 1302 in response to commands sent to the actuators. When commanded, the actuators apply forces to spacecraft 1302 that increase or decrease the velocity of spacecraft 1302 and thus translate its position, and the actuators also apply torques to spacecraft 1302 when commanded, which cause spacecraft 1302 to rotate, thereby changing its attitude or orientation. As used herein, the operation of spacecraft 1302 is determined by the operation of the actuators, which determine the motion of spacecraft 1302 that changes these quantities.

[0226] Spacecraft 1302 travels in space around, between, or near one or more gravitational bodies (such as Earth 1361, the Moon, and / or other celestial bodies, planets, stars, asteroids, comets) along an open or closed orbital path 1360. Typically, a desired or target position 1365 is given along the orbital path. A reference frame 1370 is attached to the desired position, where the origin of the frame (i.e., the zero coordinates in that reference frame) is always the coordinate of the desired position.

[0227] Spacecraft 1302 is subjected to various disturbance forces 1314. These disturbance forces may include forces not considered when determining the orbital path of spacecraft 1302. These disturbance forces act on spacecraft 1302 to move it away from a desired position in its orbit. These forces may include, but are not limited to, gravity, radiation pressure, atmospheric drag, non-spherical central bodies, and leaking propellant. Therefore, spacecraft 1302 may move away from its target position at a certain distance 1367.

[0228] Due to disturbances, spacecraft 1302 cannot always be kept in its desired position along its orbit. Therefore, it is desirable for spacecraft 1302 to remain within a window 1366 of a specified size 364 around the desired position. To this end, spacecraft 1302 is controlled to move along any path 1380 contained within the desired target window. In this example, window 1366 has a rectangular shape, but the shape of the window can vary for different implementations.

[0229] Spacecraft 1302 also frequently needs to maintain a desired orientation. For example, the spacecraft's fixed reference frame 1374 needs to be aligned with a desired reference frame, such as an inertial reference frame 1371 fixed relative to a distant star 1372, or a reference frame 1373 always oriented with a point pointing towards Earth. However, depending on the shape of spacecraft 1302, different disturbance forces 1314 can act non-uniformly on spacecraft 1302, generating disturbance torques that cause spacecraft 1302 to rotate away from its desired orientation. To compensate for the disturbance torques, momentum exchange devices 1351, such as reaction wheels, are used to absorb the disturbance torques, thereby allowing the spacecraft to maintain its desired orientation.

[0230] To prevent the momentum exchange devices from saturating and losing their ability to compensate for disturbance torque, the momentum they store must be unloaded, for example, by reducing the rotational speed of the reaction wheel. Unloading the momentum exchange devices introduces unwanted torque into spacecraft 1302. This unwanted torque is also compensated for by the thrusters.

[0231] In some implementations, the predictive controller determines the inputs to spacecraft 1302 based on a mixed integer control solution, wherein the inputs to spacecraft 1302 drive one or a combination of thrusters and momentum exchange devices, and discrete optimization variables are used to model one or a combination of discrete control decisions, switching in system dynamics, integer values ​​of thruster commands, and obstacle avoidance constraints.

[0232] In some implementations, spacecraft 1302 may be modeled as a hybrid system and a predictive controller, such as a mixed-integer model predictive controller, may be used to compute commands sent to the actuators. For example, in some implementations, commands sent to thruster 1350 may only take a discrete set of values, thus resulting in a set of binary or integer control input variables for each stage within the mixed-integer control process.

[0233] In some implementations, the predictive controller is designed to keep spacecraft 1302 outside a specific region 1385 of a specified size, approaching the desired position along its orbit. This latter region can be time-fixed or time-varying, and is often referred to as the exclusion region 1385. For this purpose, the corresponding logical inequality constraints can be modeled using an additional set of binary or integer control input variables for each stage within the mixed-integer control program. In this example, the exclusion region 1385 has a rectangular shape and is located in the corner of the desired window 1366; however, the shape and position of the exclusion region within the desired target window can vary for different implementations.

[0234] Figure 14A A schematic diagram of a vapor compression system 1400 controlled by a controller 1460 according to some embodiments is illustrated. The controller 1460 includes a predictive controller, such as a controller implementing model predictive control (MPC). Components of the vapor compression system (VCS) 1400 may include an indoor heat exchanger 1420 located in an indoor space or area 1450, an outdoor unit heat exchanger 1430 located in the surrounding environment, a compressor 1410, and an expansion valve 1440. A heat load 1415 acts on the indoor space or area 1450.

[0235] Additionally, the VCS 1400 may include a reverse-flow valve 1455 for directing high-pressure refrigerant leaving the compressor to the outdoor or indoor unit heat exchanger, and low-pressure refrigerant returning from the indoor or outdoor unit heat exchanger to the compressor inlet. When high-pressure refrigerant is directed to the outdoor unit heat exchanger, the outdoor unit heat exchanger acts as a condenser and the indoor unit heat exchanger acts as an evaporator, where the system dissipates heat from that area to the surrounding environment; this is operationally referred to as "cooling mode." Conversely, when high-pressure refrigerant is directed to the indoor unit heat exchanger, the indoor unit heat exchanger acts as a condenser and the outdoor unit heat exchanger acts as an evaporator, extracting heat from the surrounding environment and pumping that heat into the area; this is operationally referred to as "heating mode."

[0236] Figure 14B An example configuration of signals, sensors, and controllers used in a VCS 1400 according to some embodiments is illustrated. Controller 1460 reads information from sensor 1470, which is configured to measure various temperatures, pressures, flow rates, or other information about the system's operation, including measurable disturbances such as ambient air temperature. Controller 1460 may be provided with a setpoint 1466, representing the desired value of a measured signal of the process, such as the desired zone temperature. The setpoint information may come from a thermostat, a wireless remote control, or internal memory or storage medium. The controller then calculates control inputs such that some measured outputs are driven to their setpoints. These control inputs may include indoor unit fan speed 1480, outdoor unit fan speed 1481, compressor rotation speed 1482, expansion valve position 1483, and flow reversing valve position 1484. In this way, the controller controls the operation of the vapor compression system to achieve the setpoint value in the presence of disturbances 1468 acting on the system (such as heat load).

[0237] In some implementations, the VCS 1400 can be modeled as a hybrid system, and a predictive controller (such as a mixed-integer model predictive controller) is used to compute commands sent to the actuators. For example, in some implementations, commands sent to valves and / or fans may only take a discrete set of values, thus resulting in a set of binary or integer control input variables for each stage within the mixed-integer control program.

[0238] In some implementations, the predictive controller determines the inputs to the vapor compression system based on a mixed integer control solution, wherein the inputs to the vapor compression system include one or a combination of indoor unit fan speed, outdoor unit fan speed, compressor rotational speed, expansion valve position, and flow reversal valve position, and discrete optimization variables are used to model one or a combination of integer values ​​of discrete control decisions, switching in system dynamics, and commands sent to valves and / or fans.

[0239] In some implementations, the dynamic behavior of the VCS 1400 can change rapidly or even switch at certain times based on the current state of the system and the current control input values. The resulting hybrid VCS 1400 with switching dynamics can be modeled using an additional set of binary or integer control input variables for each stage within the hybrid integer control program.

[0240] Figure 15 A method 1500 for a control system according to an example embodiment is illustrated. In step 1501, the method includes receiving a feedback signal including a measurement result of the system's state. In step 1503, the method includes solving a mixed-integer optimal control optimization problem using branch-and-bound (B&B) optimization that searches for a global optimum within a search space to generate a control signal, wherein the B&B optimization iteratively partitions the search space into a nested tree of regions, and at least one region is removed from the nested tree before finding a local optimum for each region if the dual objective value of the projection of the suboptimal dual solution estimate of each region into the dual feasible space is greater than an upper bound or less than a lower bound of the global optimum maintained by the B&B optimization. In step 1505, the system is controlled based on the control signal to change the state of the system.

[0241] The embodiments described above can be implemented in any of a variety of ways. For example, the embodiments can be implemented using hardware, software, or a combination thereof. When implemented in software, the software code can execute on any suitable processor or on a set of processors, whether located in a single computer or distributed across multiple computers. Such a processor can be implemented as an integrated circuit having one or more processors within an integrated circuit assembly. However, the processor can be implemented using circuitry of any suitable format.

[0242] Furthermore, the various methods or processes outlined herein can be encoded as software that executes on one or more processors employing any of a variety of operating systems or platforms. Additionally, this software can be written using any of a large number of suitable programming languages ​​and / or programming or scripting tools, and can also be compiled into executable machine language code or intermediate code that executes on an architecture or virtual machine. Typically, the functionality of program modules can be combined or distributed as desired in various implementations.

[0243] Furthermore, embodiments of the present invention can be embodied as a method, examples of which have been provided. Actions performed as part of this method can be ordered in any suitable manner. Therefore, embodiments can be constructed in which actions are performed in a different order than those shown, which may include performing some actions simultaneously, even if these actions are shown as sequential actions in the exemplary embodiments.

[0244] Although the invention has been described by way of example of preferred embodiments, it should be understood that various other adaptations and modifications can be made within the spirit and scope of the invention. Therefore, the appended claims are intended to cover all such variations and modifications that fall within the true spirit and scope of the invention.

Claims

1. A predictive controller for a control system, the predictive controller comprising: At least one processor; and a memory having instructions stored thereon, which, when executed by the at least one processor, cause the prediction controller to: Receive feedback signals including measurement results of the system's state; A mixed-integer optimal control optimization problem is solved using branch-and-bound (B&B) optimization, which searches for a global optimum within a search space, to generate a control signal. The B&B optimization iteratively partitions the search space into a nested tree of regions, and removes at least one region from the nested tree before finding a local optimum if the dual objective value of the projection of the suboptimal dual solution estimate of each region into the dual feasible space is greater than an upper bound or less than a lower bound of the global optimum maintained by the B&B optimization. The system is controlled based on the control signal to change the state of the system. The nested region tree is formed by different convex relaxations of the integer variables in the mixed integer optimal control optimization problem. Furthermore, the local optimum within each region of the nested region tree is searched through multiple iterations using an interior-point method, such that intermediate iterations of the interior-point method generate estimates of the suboptimal dual solution for each region. The local optimum is based on an early termination process, which allows the interior point method to terminate its iterative process before finding a local optimum in the at least one region if: After one or more iterations of the interior-point method, a projection step is performed to compute estimates of dual possible solutions, and the corresponding dual objective value is greater than the upper bound of the global optimum maintained by the B&B optimization. The projection step involves solving an equality-constrained optimization problem to update the estimate of the suboptimal dual solution. This equality-constrained optimization problem enforces dual feasibility constraints and minimizes the norm of the updated Lagrange multiplier variables. The projection step applies one or more additional optimality conditions based on the positive definite weighting matrix in the objective function, and minimizes the weighted norm of the updated optimization variables.

2. The predictive controller according to claim 1, wherein, The early termination procedure also allows the interior point method to terminate its iterative process before finding a local optimum for the at least one region if: The dual feasible solution estimate is computed after one or more iterations of the interior point method, and the dual objective value is greater than the upper bound of the global optimum maintained by the B&B optimization. or The interior point method detects, after one or more iterations, that no feasible solution exists for the at least one region.

3. The predictive controller according to claim 1, wherein, The processor is configured to: When the dual objective value estimated by the suboptimal dual solution is greater than the upper limit, the projection of the suboptimal dual solution estimate used in the early termination process is performed; or The iterative process continues from the estimation of the suboptimal dual solution to find the global optimum. The local optimal solution of the at least one region is used to initialize the iterative process in one or more other regions within the nested region tree used for the B&B optimization.

4. The predictive controller according to claim 1, wherein, The projection step calculates the suboptimal dual solution estimate from the local optimum to perform the early termination process for one or more other regions within the nested region tree used for the B&B optimization.

5. The predictive controller according to claim 1, wherein, The positive definite weighting matrix in the objective of the projection step corresponds to the weighted Hessian matrix as an update of the optimization variables, the diagonal matrix with small positive weights for updating the Lagrange multipliers corresponding to the equality constraints in the convex relaxation, and the positive weights corresponding to the ratio between the relaxation variables and the Lagrange multiplier values ​​corresponding to each inequality constraint after one or more iterations in the convex solution process.

6. The predictive controller according to claim 1, wherein, The projection step updates the estimate of the suboptimal dual solution by solving the block-structured linear system, and the solution of the block-structured linear system forms the optimal solution to the equality-constrained optimization problem. Specifically, the block-structured linear system is solved by solving the simplified linear system to calculate the updates to the optimization variables, and then solving the diagonal linear system to calculate the updates to the Lagrange multiplier variables. Here, given the Chollisky decomposition from a prior projection step or from a prior iteration of the convex solution process, the Chollisky decomposition is computed or updated to compute the solution of the simplified linear system with a positive definite KKT matrix, and The KKT matrix of the simplified linear system has a block tridiagonal sparse structure due to the block-structured sparsity of the mixed integer optimal control optimization problem, and the block tridiagonal Choreski decomposition is computed or updated to compute the solution of the simplified linear system with a positive definite block tridiagonal KKT matrix.

7. The predictive controller according to claim 1, wherein, The predictive controller is implemented using Mixed Integer Model Predictive Control (MI-MPC), wherein the MI-MPC calculates the control signal based on the current state of the system and the control command, and wherein the MI-MPC calculates the control solution by solving a constrained mixed integer optimization problem at each control time step, the control solution comprising a sequence of future optimal discrete and continuous control inputs over the prediction time history of the system.

8. The predictive controller according to claim 1, wherein, The system is a vehicle, and the predictive controller determines the inputs to the vehicle based on a mixed-integer control solution, wherein the inputs to the vehicle include one or a combination of the vehicle's acceleration, the vehicle's engine torque, braking torque, and steering angle, and discrete optimization variables are used to model one or a combination of discrete control decisions, switching, gear shifting, and obstacle avoidance constraints in system dynamics.

9. The predictive controller according to claim 1, wherein, The system is a spacecraft, and the predictive controller determines the input to the spacecraft based on a mixed integer control solution, wherein the input to the spacecraft drives one or a combination of thrusters and momentum exchange devices, and discrete optimization variables are used to model one or a combination of discrete control decisions, switching in system dynamics, integer values ​​of thruster commands, and obstacle avoidance constraints.

10. The predictive controller according to claim 1, wherein, The system is a vapor compression system, and the predictive controller determines the inputs to the vapor compression system based on a mixed integer control solution, wherein the inputs to the vapor compression system include one or a combination of indoor unit fan speed, outdoor unit fan speed, compressor rotational speed, expansion valve position, and flow reversal valve position, and discrete optimization variables are used to model one or a combination of discrete control decisions, switching in system dynamics, and integer values ​​of commands sent to valves and / or to fans.

11. A method for a control system, wherein, This method uses: At least one processor; as well as A memory having instructions stored thereon, the instructions implementing the steps of the method when executed by the at least one processor, the steps including: Receive feedback signals including measurement results of the system's state; A mixed-integer optimal control optimization problem is solved using branch-and-bound (B&B) optimization, which searches for a global optimum within the search space, to generate a control signal. The B&B optimization iteratively partitions the search space into a nested tree of regions, and removes at least one region from the nested tree before finding a local optimum if the dual objective value of the projection of the suboptimal dual solution estimate of each region into the dual feasible space is greater than an upper bound or less than a lower bound of the global optimum maintained by the B&B optimization. The system is controlled based on the control signal to change the state of the system. The nested region tree is formed by different convex relaxations of the integer variables in the mixed integer optimal control optimization problem. Furthermore, the local optimum within each region of the nested region tree is searched through multiple iterations using an interior-point method, such that intermediate iterations of the interior-point method generate estimates of the suboptimal dual solution for each region. The local optimum is based on an early termination process, which allows the interior point method to terminate its iterative process before finding a local optimum in the at least one region if: After one or more iterations of the interior-point method, a projection step is performed to compute a dual feasible solution estimate, and the corresponding dual objective value is greater than the upper bound of the global optimum maintained by the B&B optimization. The projection step involves solving an equality-constrained optimization problem to update the estimate of the suboptimal dual solution. This equality-constrained optimization problem enforces dual feasibility constraints and minimizes the norm of the updated Lagrange multiplier variables. The projection step applies one or more additional optimality conditions based on the positive definite weighting matrix in the objective function, and minimizes the weighted norm of the updated optimization variables.

12. The method of claim 11, wherein, The early termination procedure also allows the interior point method to terminate its iterative process before finding a local optimum for the at least one region if: After one or more iterations of the interior point method, a dual feasible solution estimate is computed, and the dual objective value is greater than the upper bound of the global optimum maintained by the B&B optimization. or The interior point method detects, after one or more iterations, that no feasible solution exists for the at least one region.

13. A non-transitory computer-readable storage medium on which a program executable by a processor for performing a method comprising the following operations is implemented: It accepts feedback signals, including measurement results of the system's state; Branch and bound (B&B) optimization, which searches for the global optimum within the search space, is used to solve the mixed-integer optimal control optimization problem to generate control signals, where... The B&B optimization iteratively partitions the search space into a nested tree of regions, and when the dual objective value of the projection of the suboptimal dual solution estimate of each region into the dual feasible space is greater than the upper limit of the global optimum maintained by the B&B optimization or less than the lower limit of the global optimum, at least one region is removed from the nested tree of regions before finding the local optimum of each region. as well as The system is controlled based on the control signal to change the state of the system. The nested region tree is formed by different convex relaxations of the integer variables in the mixed integer optimal control optimization problem. Furthermore, the local optimum within each region of the nested region tree is searched through multiple iterations using an interior-point method, such that intermediate iterations of the interior-point method generate estimates of the suboptimal dual solution for each region. The local optimum is based on an early termination process, which allows the interior point method to terminate its iterative process before finding a local optimum in the at least one region if: After one or more iterations of the interior-point method, a projection step is performed to compute a dual feasible solution estimate, and the corresponding dual objective value is greater than the upper bound of the global optimum maintained by the B&B optimization. The projection step involves solving an equality-constrained optimization problem to update the estimate of the suboptimal dual solution. This equality-constrained optimization problem enforces dual feasibility constraints and minimizes the norm of the updated Lagrange multiplier variables. The projection step applies one or more additional optimality conditions based on the positive definite weighting matrix in the objective function, and minimizes the weighted norm of the updated optimization variables.

Images