Sequential Convexification Method for Model Predictive Control of Nonlinear Systems with Continuous and Discrete Elements of Operations
To control a hybrid dynamical system, a predictive feedback controller formulates a mixedinteger nonlinear programming (MINLP) problem including nonlinear functions of continuous optimization variables representing the continuous elements of the operation of the hybrid dynamical system and/or one or multiple linear functions of integer optimization variables representing the discrete elements of the operation of the hybrid dynamical system. The MINLP problem is formulated into a separable format ensuring that the discrete elements of the operation are present only in the linear functions of the MINLP problem. The MINLP problem is solved over multiple iterations using a partial convexification of a portion of a space of the solution including a current solution guess. The partial convexification produces a convex approximation of the nonlinear functions of the MINLP without approximating the linear functions of the MINLP to produce a partially convexified MINLP.
Latest Mitsubishi Electric Research Laboratories, Inc. Patents:
 System and Method for Robust Pivoting for ReOrienting Parts during Robotic Assembly
 Method and System for Fault Location and Protection of InverterDominated Islanded Ungrounded Microgrids
 ModelBased Control with Uncertain Motion Model
 Motor Assembly for Linear DirectDrive Motor
 System and Method for Detecting BrokenBar Fault in SquirrelCage Induction Motors
The present disclosure relates generally to mixedinteger nonlinear optimizationbased control, and more particularly to a sequential convexification method and apparatus for model predictive control of systems that are described by nonlinear dynamics with continuous and discrete elements of operations.
BACKGROUNDOptimization based decision making, planning and control techniques, such as model predictive control (MPC), allow a modelbased design framework in which the system dynamics, the system requirements and constraints can directly be taken into account. This framework has been extended to hybrid dynamical systems, including both continuous and discrete decision variables, which provides a powerful technique to model a large range of problems, e.g., including dynamical systems with mode switchings or systems with quantized actuation, problems with logic rules, temporal logic specifications or obstacle avoidance constraints. However, the resulting optimization problems are highly nonconvex, and therefore difficult to solve in practice, because they contain variables which only take discrete values (e.g., binary or integer values). When using nonlinear system dynamics, one or multiple nonlinear constraint functions and/or a nonlinear objective function, the resulting optimal control problem (OCP) can be formulated as a mixedinteger nonlinear programming (MINLP) problem, which is NPhard and therefore computationally difficult to solve.
In general, mixedinteger nonlinear model predictive control (MINMPC) requires the solution of a nonconvex MINLP, i.e., the optimization problem is nonconvex even after relaxing the integrality constraints, e.g., due to the nonlinear system dynamics, nonlinear constraint functions and/or objective functions. Most successful global optimization algorithms for MINLPs require convexity of the objective and constraint functions, which can therefore not be used to solve MINMPC problems to global optimality. Even though global optimization algorithms exist for nonconvex MINLPs, e.g., using relaxations of factorable problems, they are usually computationally very expensive and hence generally not yet practical for realtime implementations of MINMPC.
Decision making, planning or control for hybrid systems aims to solve an MINLP at every sampling time instant to enable realtime MINMPC applications. Hence, some methods therefore focus on approximate or heuristic techniques to find feasible but (possibly) suboptimal solutions to MINLPs within strict timing requirements. Some existing techniques are based on global algorithms for convex MINLPs, which can be used to find approximate solutions to nonconvex MINLPs, e.g., using outer approximation or hybrid branchandbound (hB&B) methods. Specifically for nonconvex MINMPC, a variant of the realtime iterations (RTI) algorithm has been proposed based on outer convexification in combination with rounding schemes. However, when inequality constraints depend directly on the discrete decision variables, the latter approach requires solving mathematical programs with vanishing constraints, which are particularly challenging.
Sequential convexification techniques, e.g., using sequential convex programming (SCP) or sequential quadratic programming (SQP) methods form a popular technique to solve general nonlinear programming (NLP) problems. In particular, sequential convexification techniques have been used successfully for the realtime implementation of nonlinear model predictive control (NMPC) with smooth nonlinear dynamics, nonlinear constraint functions and/or inequality constraints. However, there is a need for the extension of these methods to NMPC for systems with both continuous and discrete elements of operation, i.e., including continuous as well as integer and/or binary decision variables in MINMPC.
In recent prior work, a mixedinteger sequential quadratic programming (MISQP) method was proposed that is based on the use of a trust region radius for both continuous and integer optimization variables. The method requires the solution of a mixedinteger quadratic programming (MIQP) subproblem at each iteration, which can be solved efficiently, e.g., using state of the art branchandbound (B&B) optimization methods. The standard MISQP method however relies on the assumption that integer variables have a smooth influence on the MINLP, i.e., incrementing an integer variable by one leads to a small change of function values. However, the latter assumption is generally not true for MINMPC because, for example, the constrained optimization problem may include binary variables that have a large influence on the optimal control trajectories.
Therefore, there is a need for a sequential convexification method that is more generally applicable to MINMPC problems, which is the aim of the system and method described in the present invention.
SUMMARYEmbodiments of the invention are based on the solution of a sequence of one or multiple mixedinteger convex programming (MICP) subproblems, in which the preparation of each subproblem is performed based on a partial convexification technique in order to compute a feasible but possibly suboptimal solution of the mixedinteger nonlinear optimal control problem at each sampling time instant of the proposed MINMPC controller. The solution of each MICP subproblem can be used to compute an update to the optimal solution guess for all integer and/or binary decision variables, as well as a new search direction for the continuous decision variables. In addition, based on the updated values for integer and/or binary decision variables, and based on the new search direction for continuous decision variables, the current solution guess for the continuous decision variables can be updated in each iteration of the mixedinteger sequential convex programming (MISCP) optimization algorithm.
Some embodiments of the invention are based on the realization that any MINLP can be reformulated as a different but mathematically equivalent MINLP in separable format, in which all integer and/or binary decision variables enter linearly in all constraint and objective functions. The latter reformulation can be achieved, for example, by defining one or multiple auxiliary continuous optimization variables to ensure that all integer and/or binary optimization variables enter linearly in constraint and objective functions. Specifically, all nonlinear functions that may be present in the constraint and objective functions of the MINMPC formulation in separable format therefore depend only on continuous optimization variables. Some embodiments of the invention are based on the realization that the latter linear dependency of the constraint and objective functions in the MINMPC formulation on all integer and/or binary decision variables can be used to avoid the smoothness requirement, i.e., incrementing an integer variable by one leads to a small change of function values, which limits the applicability of standard MISQP algorithms for MINMPC.
A partial convexification technique is used to prepare the MICP subproblem in each iteration of the proposed MISCP optimization algorithm. Some embodiments of the invention are based on the realization that the partial convexification technique needs to be applied only to smooth nonlinear functions, due to the linear dependency of the constraint and objective functions in the MINMPC formulation on all integer and/or binary decision variables. In some embodiments of the invention, the partial convexification technique is based on a local linearization of all smooth nonlinear functions that may be present in the constraint and objective functions of the MINMPC formulation, based on a solution guess for the continuous optimization variables in the current iteration of the MISCP algorithm. Alternatively, in other embodiments of the invention, a partial convexification technique can be used to compute more general convex approximations of nonlinear functions in one or multiple inequality constraints of the MINMPC problem formulation, resulting in convex quadratic inequality constraints in a mixedinteger quadratically constrained quadratic programming (MIQCQP) subproblem, or convex second order cone constraints in a mixedinteger second order cone programming (MISOCP) subproblem.
Embodiments of the invention are based on the realization that the solution of an MICP subproblem in each iteration of the proposed MISCP algorithm can be computed relatively fast, thanks to progress that has been made over past decades in the development of state of the art solvers for MICPs. For example, branchandbound methods can be used to efficiently solve a mixedinteger quadratic programming (MIQP) or a mixedinteger linear programming (MILP) subproblem. State of the art branchandbound methods for MIQPs and/or MILPs may include advanced primal heuristics, branching strategies, presolve operations, cut generation techniques, and convex solvers with early termination and infeasibility detection to effectively reduce the size of the branchandbound search tree and therefore reduce the amount of convex relaxations that need to be solved in order to compute the globally optimal solution to each MIQP/MILP.
Some embodiments of the invention are based on the realization that the update of the continuous decision variables in each MISCP iteration can be performed in a particular way to ensure some amount of progress is made in computing a feasible but possibly suboptimal solution to the MINLP problem. Some embodiments of the invention use a globalization strategy based on a merit function that quantifies a combination of optimality and constraint satisfaction for a solution guess of values for continuous and discrete decision variables. One example of a merit function is based on an l_{1 }penalty function, applied to each of the equality constraints and to the violation of each of the inequality constraints in the MINMPC problem formulation, excluding the integrality constraints which are automatically satisfied by design for each of the integer and/or binary decision variables at each iteration of the proposed MISCP optimization algorithm.
In some embodiments of the invention, the globalization strategy is based on a line search method that computes a step size α^{k}∈(0,1] in order to update the continuous optimization variables y^{k+1}=y^{k}+α^{k}Δy^{k}, based on the search direction Δy^{k }for the continuous optimization variables y that is computed in the kth iteration of the MISCP optimization algorithm. The step size selection aims to satisfy a condition for the sufficient decrease of the value for the merit function, evaluated at the new solution guess y^{k+1}, compared to the original value of the merit function for the solution guess y^{k }and given a directional derivative of the merit function at each iteration of the MISCP algorithm.
In some embodiments of the invention, the globalization strategy is based on a trustregion method that computes a subregion in the space of continuous decision variables, in which the local convex approximation of the nonlinear constraint and/or objective functions is sufficiently accurate. The accuracy of the partial convexification can be approximately evaluated based on the ratio of the actual versus the predicted reduction in the value of the merit function from one solution guess to the next in each iteration of the MISCP optimization algorithm. Some embodiments of the invention are based on a trustregion radius, whose value can be either increased, decreased or kept the same in each iteration of the algorithm.
Some embodiments of the invention are based on the realization that a relatively small trustregion radius may be used when the MISCP optimization algorithm is sufficiently close to a feasible but possibly suboptimal solution to the MINMPC problem, resulting in a considerable reduction of the computational cost to solve the MICP subproblem due to the relatively small trustregion radius.
Some embodiments of the invention are based on a warm starting strategy that computes a guess for the optimal values of the continuous and of the integer and/or binary decision variables, based on the feasible but possibly suboptimal solution to the MINMPC problem at the previous sampling time instant. For example, in some embodiments of the invention, a time shifting procedure by one sampling time period can be used to warm start the solution guess for the MISCP optimization algorithm, given the optimal or approximately optimal solution to the MINLP at the previous sampling time instant of the proposed MINMPC controller.
Some embodiments of the invention are based on a limited number of one or multiple solutions of MICP subproblems, after which some or all of the integer and/or binary decision variables remain fixed, resulting in a number of one or multiple solutions of convex programming (CP) subproblems. The upper bound on the number of MICP solutions can be chosen sufficiently large to enable the MISCP optimization algorithm to compute a good feasible but possibly suboptimal solution to the MINMPC problem, while remaining relatively small such that the computational cost of the algorithm can be reduced considerably. For example, in some embodiments of the invention, only a single MICP subproblem needs to be solved at each sampling time instant of the realtime feasible MINMPC controller.
Some embodiments of the invention are based on the realization that one or multiple solutions of CP subproblems can be performed with a computational cost that is considerably smaller than one or multiple solutions of MICP subproblems. Some embodiments of the invention are based on the realization that fixing of some or all of the integer and/or binary decision variables, after a limited number of one or multiple solutions of MICP subproblems, may prevent cycling of the MISCP optimization algorithm in the proposed MINMPC controller.
In some embodiments of the invention, the number of MISCP iterations is determined by a termination condition of the optimization algorithm, which can be based on a norm of the KarushKuhnTucker (KKT) necessary conditions of optimality for the MINLP, excluding the integrality conditions.
In some embodiments of the invention, a homotopytype penalty method can be used to adjust the cost function in each MICP subproblem to increasingly enforce the MISCP algorithm to compute an update to the optimal solution guess for some or all of the integer and/or binary decision variables that remains close to the solution guess of the integer and/or binary decision variables at the previous MISCP iteration. For example, a term w_{i}z_{i }can be added to the cost function of the MICP minimization subproblem, where w_{i}>0 is a positive weight value, to ensure a binary variable z_{i}∈{0,1} to remain close to a solution guess z_{i}^{k}=0. Alternatively, a term w_{i}(1−z_{i}) can be added to the cost function of the MICP minimization subproblem, where w_{i}>0 is a positive weight value, to ensure a binary variable z_{i}∈{0,1} to remain close to a solution guess z_{i}^{k}=1.
Some embodiments of the invention are based on the realization that the use of a homotopytype penalty method for some or all of the integer and/or binary decision variables may prevent cycling in the MISCP optimization algorithm. In addition, some embodiments of the invention are based on the realization that the use of a homotopytype penalty method for some or all of the integer and/or binary decision variables may considerably reduce the computational cost of solving the MICP subproblems in the MISCP optimization algorithm.
Accordingly, one embodiment discloses a predictive feedback controller for controlling a hybrid dynamical system with nonlinear dynamics and continuous and discrete elements of operation, the predictive feedback controller comprising: at least one processor; and a memory having instructions stored thereon that, when executed by the at least one processor, cause the predictive feedback controller to:

 accept feedback signal including measurements indicative of a current state of the hybrid dynamical system including one or a combination of a current state of the predictive controller, a current state of one or multiple actuators of the hybrid dynamical system, and a current state of outputs of the hybrid dynamical system;
 formulate a mixedinteger nonlinear programming (MINLP) problem optimizing an objective function subject to one or multiple constraints with a solution indicative of a control command for changing the current state of the hybrid dynamical system according to a control objective, wherein the constraints include equality constraints, inequality constraints, or both, and wherein the constraints and the control objective of the MINLP problem include one or multiple nonlinear functions of continuous optimization variables representing the continuous elements of the operation of the hybrid dynamical system and one or multiple linear functions of integer optimization variables representing the discrete elements of the operation of the hybrid dynamical system, such that the MINLP problem is formulated into a separable format ensuring that the discrete elements of the operation are present only in the linear functions of the MINLP problem;
 solve the MINLP problem over multiple iterations of a sequential convexificationbased optimization procedure until a termination condition is met, wherein, to perform an iteration, the predictive feedback controller is configured to perform a partial convexification of a portion of a space of the solution including a current solution guess, wherein the partial convexification produces a convex approximation of the nonlinear functions of the MINLP without approximating the linear functions of the MINLP to produce a partially convexified MINLP; and update the current solution guess by solving a mixedinteger convex programming (MICP) formulation of the partially convexified MINLP problem; and
 submit the control command generated according to the solution of the MINLP problem to the hybrid dynamical system thereby causing a change of the current state of the hybrid dynamical system.
Some embodiments of the present disclosure provide a system and a method for controlling an operation of a system or a system using a predictive controller. An example of the predictive controller is a mixedinteger nonlinear model predictive controller (MINMPC) determining control inputs based on a nonlinear model of the controlled system having continuous and discrete elements of operations.
The system 120, as referred herein, can be any machine or device controlled by certain manipulation input signals, e.g., control signal 111 (inputs). The control input signal can possibly include continuous elements such as voltages, pressures, forces, torques, steering angles, velocities, and temperatures, as well as discrete elements such as energy levels, quantized valve inputs, gear shifts, on/off actuation, lane selection, and obstacle avoidance decision variables. The system 120 returns some controlled output signals 103 (outputs), possibly including continuous elements such as currents, flows, velocities, positions, temperatures, heading and steering angles, as well as discrete elements such as energy levels, quantized valve states, gear status, on/off status, and lane position, etc. The output values are related in part to previous output values of the system, and in part to previous and current input values. The dependency on previous inputs and previous outputs is encoded in the state of the system. The operation of the system, e.g., a motion of components of the system, can include a sequence of output values generated by the system following the application of certain input values.
The system model 102 may include a set of mathematical equations that describe how the system outputs change over time as functions of current and previous inputs, and the previous outputs. The mathematical equations can include one or multiple smooth equations depending on continuous variables as well as one or multiple mixedinteger equations depending on both continuous and discrete variables. Each function in the mathematical equations can be either a linear or a smooth nonlinear function. The state of the system 120 is any set of information, in general time varying, for instance an appropriate subset of current and previous inputs and outputs, that, together with the model of the system and future inputs, can uniquely define the future motion of the system.
The system 120 can be subject to physical limitations and specification constraints 104 limiting the range where the outputs, the inputs, and also possibly the states of the system 120 are allowed to operate. The constraints can include one or multiple smooth equations depending on continuous variables as well as one or multiple mixedinteger equations depending on both continuous and discrete variables. The constraint functions can be either linear or smooth nonlinear functions. Some embodiments of the invention are based on the realization that first or higher order directional derivatives can be computed for each of the smooth nonlinear functions in the constraints 104.
The predictive controller 110 can be implemented in hardware or as a software program executed in a processor, e.g., a microprocessor, which at fixed or variable control period sampling intervals receives the estimated state 121 of the system 120 and the desired motion command 101 and determines, using this information, the inputs, e.g., the control signal 111, for operating the system 120. The predictive controller 110 further solves an optimal control structured mixedinteger nonlinear programming (MINLP) problem using a mixedinteger sequential convex programming (MISCP) method, optimizing a current solution of the objective function subject to one or multiple equality and/or inequality constraints over multiple iterations until a termination condition is met, e.g., until the solution is feasible and (locally) optimal for the optimal control structured MINLP. Each iteration of the MISCP method performs a partial convexification of a portion of a space of the solution including the current solution, wherein the partial convexification produces a convex approximation of the smooth nonlinear functions of the MINLP without approximating the linear functions of the MINLP to produce a partially convexified MINLP. Then, each iteration can update the current solution by solving a mixedinteger convex programming (MICP) formulation of the partially convexified MINLP problem to produce a control signal 111. The predictive controller 110, further, controls the system 120 based on the control signal 111 to change the state of the system 120.
The estimator 130 can be implemented in hardware or as a software program executed in a processor, either the same or a different processor from the controller 110, which at fixed or variable control period sampling intervals receives the outputs of the system 103 and determines, using the new and the previous output measurements, the estimated state 121 of the system 120.
Thus, by using the MISCP optimization method, the processor computes a feasible and (locally) optimal solution for the optimal control structured MINLP by solving a sequence of MICP subproblems, where the formulation of each MICP subproblem is based on a partial convexification of a portion of the solution space. Due to the reduced computational complexity of solving the MICP subproblems compared to the computational complexity of solving the original MINLP, the processor is enabled to accurately determine feasible and (locally) optimal solutions to control the state of the system 120 in less time. Accordingly, the processor achieves fast processing speed with high accuracy.
An example of a tracking controller 115 can be based on a proportionalintegralderivative (PID) controller to track the timevarying reference motion trajectory that is computed by the MINMPC controller 110. In some embodiments of the invention, the tracking controller 115 can be based on a model predictive controller (MPC) with a relatively simplified dynamical model and a set of simplified constraints, and therefore requiring a relatively small computational complexity. For example, the MPC tracking controller 115 can be based on a linearquadratic objective function, linear dynamical model and linear equality and inequality constraints, which requires the solution of a convex linear programming (LP) or a convex quadratic programming (QP) problem that is computationally much easier to solve than the MINLP that is solved by the MINMPC controller 110.
Some embodiments of the invention are based on the realization that a relatively long prediction horizon can be used for the mixedinteger nonlinear predictive controller 110 and a relatively short prediction horizon can be used for the tracking controller 115, such that a relatively fast sampling rate can be used for the tracking controller, e.g., with a sampling period of 10100 milliseconds, while a relatively slow sampling rate can be used for the mixedinteger nonlinear predictive controller 110, e.g., with a sampling period of 0.52 seconds.
Some embodiments of the invention are based on the realization that a relatively highlevel, lowaccuracy dynamical model 102 and constraints 104 can be used in the MINLP formulation of the mixedinteger nonlinear predictive controller 110, due to the computational complexity of solving MINLPs, while a relatively lowlevel, highaccuracy dynamical model 102 and constraints 104 can be used in the design and formulation of the computationally cheap tracking controller 115. [what could be the minimum ratio of long/short perdition horizons? How the partial convexification helps to achieve this ratio?]
In some embodiments of the invention, the solution of this inequality constrained MINLP problem 250 uses one or multiple state and control values over a prediction time horizon, and potentially other MINLP solution information from the previous control time step 210, which can be read from the memory. This concept is called warm or hotstarting of the optimization algorithm and it can reduce the required computational effort of the MINMPC controller in some embodiments. In a similar fashion, the corresponding solution vector 255 can be used to update and store in memory a sequence of one or multiple optimal state and control values over a prediction time horizon, and potentially other MINLP solution information for the next control time step 235.
In some embodiments, the mixedinteger optimization algorithm is based on a search algorithm to solve the MICP subproblem, which is the result of a partial convexification step in each iteration of the sequential convexification algorithm, such that the MINMPC controller updates and stores additional mixedinteger program solution information 260 in order to reduce the computational effort of the search algorithm at one or multiple iterations in the current and/or in the next control time step. In one embodiment, the MICP subproblem at each iteration is solved using a branchandbound optimization method and the warm starting information can include data related to the nodes in the binary search tree that are part of the solution path from the root node to the leaf node where the optimal integerfeasible control solution is found, in order to improve the node selection and variable branching strategies from one iteration to the next.
In some embodiments of the invention, the nonlinear dynamical model of the system 263 is described by one or multiple linear and/or smooth nonlinear differential equations. In some embodiments of the invention, the dynamical model of the system 263 describes a linear or nonlinear hybrid system with state and/or inputdependent jumps in the dynamic equations, for example, including piecewise linear and/or piecewise smooth nonlinear equations. Specifically, an equation x_{k+1}=ψ_{k}(x_{k},u_{k},w_{k}) defines the state variables at the next time step t_{k+1}, given the state variables x_{k}, the control inputs u_{k }and the integer variables w_{k }at the previous time step t_{k }within the prediction time horizon k=0, 1, . . . , N−1.
In general, the linear discrete equality constraints 266 can state that a linear function of state and control variables is constrained to be equal to one of a discrete set of values. In some embodiments of the invention, the linear discrete equality constraints 266 can include integrality constraints, for example, a constraint w_{k,j}∈ on a particular optimization variable w_{k,j }to take only integer values. In some embodiments, the linear discrete equality constraints 266 can include binary equality constraints, for example, a constraint w_{k,j}∈{0,1} on a particular optimization variable w_{k,j }to be equal to either 0 or 1.
In some embodiments of the invention, the objective function 261 can include a summation of a stage cost within the prediction time horizon k=0, 1, . . . , N−1 and a terminal cost at a final time step t_{N}. For example, in some embodiments, the stage cost l_{k}(x_{k},u_{k},w_{k}) and the terminal cost m(x_{N},w_{N}) can include linear, linearquadratic and/or nonlinear functions. In case the control command 101 includes a reference trajectory of state and control values {x_{k}^{ref},u_{k}^{ref}}_{k= . . . }, the stage and terminal cost functions could read, for example, as follows
l_{k}(x_{k},u_{k},w_{k})=∥x_{k}−x_{k}^{ref}∥_{Q}^{2}+∥u_{k}−u_{k}^{ref}∥_{R}^{2 }
m(x_{N},w_{N})=∥x_{N}−x_{N}^{ref}∥_{P}^{2 }
where the matrices Q, R, and P are typically symmetric and positive semidefinite, and ∥x_{k}−x_{k}^{ref}∥_{Q}^{2}=(x_{k}−x_{k}^{ref})^{T}Q(x_{k}−x_{k}^{ref}). Similarly, in case the control command 101 includes a reference trajectory of state, control and integer values {x_{k}^{ref},u_{k}^{ref},w_{k}^{ref}}_{k= . . . }, the stage and terminal cost in the objective function 261 could read, for example, as follows
l_{k}(x_{k},u_{k},w_{k})=∥x_{k}−x_{k}^{ref}∥_{Q}^{2}+∥u_{k}−u_{k}^{ref}∥_{R}^{2}+Σ_{j}c_{k,j}w_{k,j}−w_{k}^{ref}
m(x_{N},w_{N})=∥x_{N}−x_{N}^{ref}∥_{P}^{2}+Σ_{j}c_{N,j}w_{N,j}−w_{N}^{ref}
In some embodiments of the invention, the optimal control structured MINLP 260 can be reformulated trivially in separable format 270, for example, because the functions in the objective 261, the functions in the equality constraints 263, and the functions in the inequality constraints 264265 are defined as follows
l_{k}(x_{k},u_{k},w_{k})={tilde over (l)}_{k}(x_{k},u_{k})+c_{k}^{T}w_{k }
m(x_{N},w_{N})={tilde over (m)}(x_{N})+c_{N}^{T}w_{N }
ψ_{k}(x_{k},u_{k},w_{k})={tilde over (ψ)}_{k}(x_{k},u_{k})+D_{k}w_{k }
h_{k}(x_{k},u_{k},w_{k})={tilde over (h)}_{k}(x_{k},u_{k})+E_{k}w_{k }
h_{N}(x_{N},w_{N})={tilde over (h)}_{N}(x_{N})+E_{N}w_{N }
to define the functions in the objective 281, the functions in the equality constraints 283, and the functions in the inequality constraints 284285 in separable format, based on matrices D_{k}, E_{k}, E_{N }and vectors c_{k }for each of the time steps within the prediction time horizon k=0, 1, . . . , N−1.
In some embodiments of the invention, the optimal control structured MINLP 260 can be reformulated in separable format 270, for example, by defining one or multiple additional continuous input variables ū_{k }and/or one or multiple additional integer optimization variables
{tilde over (l)}_{k}(x_{k},ũ_{k})+c_{k}{tilde over (w)}_{k }
{tilde over (m)}(x_{N})+c_{N}^{T}{tilde over (w)}_{N }
{tilde over (ψ)}_{k}(x_{k},ũ_{k})+D_{k}{tilde over (w)}_{k }
{tilde over (h)}_{k}(x_{k},ũ_{k})+E_{k}{tilde over (w)}_{k }
{tilde over (h)}_{N}(x_{N})+E_{N}{tilde over (w)}_{N }
where
are defined. For example, in case all the integer variables w_{k }enter the system dynamics ψ_{k}(x_{k},u_{k},w_{k}) 263 nonlinearly, then the system dynamics can be defined in separable format 283 as follows
{tilde over (ψ)}_{k}(x_{k},ũ_{k})=ψ_{k}(x_{k},u_{k},ū_{k})
0=ū_{k}−w_{k }
where ū_{k,j}=w_{k,j }∀j, even though ū_{k,j}∈ and w_{k,j}∈.
Some embodiments of the invention are based on the realization that the MINLP in separable format 280 is mathematically equivalent to the original MINLP 260, i.e., the MINLP in separable format 280 is infeasible if and only if the original MINLP 260 is infeasible, and a feasible and (locally) optimal solution to the MINLP in separable format 280 can be used to construct a feasible and (locally) optimal solution of the original MINLP 260. For example, in some embodiments of the invention, an affine or nonlinear transformation exists to transform a feasible and (locally) optimal solution to the MINLP in separable format 280 into a feasible and (locally) optimal solution of the original MINLP 260.
For example, the MINLP matrices and vectors 246 can include matrices D_{k}, E_{k}, E_{N }and vectors c_{k }in the constraint and objective functions of the optimal control structured MINLP 280. Similarly, the MINLP functions 247 can include functions l_{k}(⋅), m(⋅), ψ_{k}(⋅), h_{k}(⋅) and h_{N}(⋅) in the constraint and objective functions of the optimal control structured MINLP 280.
In some embodiments of the invention, the nonlinear dynamical model of the system 283 is described by one or multiple linear and/or smooth nonlinear differential equations. In some embodiments of the invention, the dynamical model of the system 283 describes a linear or nonlinear hybrid system with state and/or inputdependent jumps in the dynamic equations, for example, including piecewise linear and/or piecewise smooth nonlinear equations. Specifically, an equation x_{k+1}=ψ_{k}(x_{k},u_{k})+D_{k}w_{k }defines the state variables at the next time step t_{k+1}, given the state variables x_{k}, the control inputs u_{k }and the integer variables w_{k }at the previous time step t_{k }within the prediction time horizon k=0, 1, . . . , N−1.
In general, the linear discrete equality constraints 286 can state that a linear function of state and control variables is constrained to be equal to one of a discrete set of values. In some embodiments of the invention, the linear discrete equality constraints 286 can include integrality constraints, for example, a constraint w_{k,j}∈ on a particular optimization variable w_{k,j }to take only integer values. In some embodiments, the linear discrete equality constraints 286 can include binary equality constraints, for example, a constraint w_{k,j}∈{0,1} on a particular optimization variable w_{k,j }to be equal to either 0 or 1.
In some embodiments of the invention, the objective function 281 can include a summation of a stage cost within the prediction time horizon k=0, 1, . . . , N−1 and a terminal cost at a final time step t_{N}. For example, in some embodiments, the stage cost l_{k}(x_{k},u_{k}) and the terminal cost m(x_{N}) can include linear, linearquadratic and/or nonlinear functions. In case the control command 101 includes a reference trajectory of state and control values {x_{k}^{ref},u_{k}^{ref}}_{k=0, . . . }, the stage and terminal cost functions could read, for example, as follows
l_{k}(x_{k},u_{k})=∥x_{k}−x_{k}^{ref}∥_{Q}^{2}+∥u_{k}−u_{k}^{ref}∥_{R}^{2 }
m(x_{N})=∥x_{N}−x_{N}^{ref}∥_{P}^{2 }
where the matrices Q, R, and P are typically symmetric and positive semidefinite, and ∥x_{k}−x_{k}^{ref}∥_{Q}^{2}=(x_{k}−x_{k}^{ref})^{T}Q(x_{k}−x_{k}^{ref}).
Each iteration of the sequential convexificationbased optimization procedure includes checking whether the current intermediate solution guess is feasible and/or (locally) optimal 350, in which case either a solution is found 255 or a new iteration is performed to update the current intermediate MINLP solution guess 335 based on a partial convexification step 315 and an MICP subproblem solution 320 and a current iteration number k can be updated 340. In some embodiments of the invention, the termination condition 350 includes checking whether a current intermediate solution guess is both feasible, i.e., the solution satisfies all equality and inequality constraints in the MINLP, whether it is sufficiently close to the globally optimal solution and/or whether the computational cost of the iterative optimization procedure has reached a particular time limit. For example, a maximum number of iterations can be imposed to ensure that the iterative optimization procedure terminates in a deterministic runtime.
Some embodiments of the invention are based on the realization that each iteration of the sequential convexificationbased optimization procedure to solve the MINLP in separable format 270 can be based on a partial convexification step 315 for only the smooth nonlinear part of the MINMPC problem in separable format, i.e., including a convexification step for each of the nonlinear equality and nonlinear inequality constraints and/or for each of the nonlinear objective functions. The partial convexification step 315 does not change the linear functions depending on one or multiple integer variables, for example, the linear constraint and objective functions in the separable MINLP format 280 based on matrices D_{k}, E_{k}, E_{N }and vectors c_{k }in the constraint and objective functions of the optimal control structured MINLP 280. Specifically, the partial convexification step 315 does not change the linear discrete equality constraints 286, resulting in a mixedinteger convex programming (MICP) subproblem to compute a search direction for continuous and integer variables (Δy^{k},Δz^{k}) 320.
Some embodiments of the invention are based on the realization that, depending on the particular implementation of the partial convexification step 315, the MICP subproblem can correspond, for example, to a mixedinteger linear programming (MILP), mixedinteger quadratic programming (MIQP), mixedinteger quadratically constrained quadratic programming (MIQCQP) or a mixedinteger second order cone programming (MISOCP) subproblem.
In some embodiments of the invention, the MICP subproblem solution (Δy^{k},Δz^{k}) 320 can be used directly to update the intermediate solution guess for integer variables z^{k+1}=z^{k}+Δz^{k }325, and it can be used to select a step size 0≤α^{k}≤1 to update the intermediate solution guess for continuous variables y^{k+1}=y^{k}+α^{k}Δy^{k }330, resulting in an updated MINLP solution guess y^{k+1}∈^{n}^{y}, z^{k+1}∈^{n}^{z }335 that can be used in the next iteration k←k+1 340. In some embodiments of the invention, a line search method can be used to select a step size 0≤α^{k}≤1 in the update of continuous variables 330, in order to ensure sufficient progress is made in each iteration of the optimization procedure with respect to one or multiple feasibility and optimality conditions for a solution to the MINLP in separable format 270. For example, a merit function can be used to quantify a combination of optimality and constraint satisfaction for a solution guess of values for continuous and discrete decision variables.
In some embodiments of the invention, the MICP subproblem solution includes (Δy^{k},Δz^{k}), subject to trustregion constraints that read as ∥MΔy∥_{p}≤d_{k}, i.e., a pnorm value of the update of continuous variables Δy^{k }is smaller than or equal to a trustregion radius value d_{k}. For example, a value of p=1, p=2, or p=∞ can be used to result in trustregion constraints in the MICP subproblem 360 based on a 1norm, 2norm, or ∞norm function, respectively. In each iteration of the sequential convexificationbased optimization procedure, the trustregion radius value d_{k }can be updated 345, based on the MICP subproblem solution 360 and the updated MINLP solution guess y^{k+1}∈^{n}^{y}, z^{k+1}∈^{n}^{z }335.
In some embodiments of the invention, the update step in each iteration of the optimization procedure is accepted or not based on one or multiple conditions which are checked 355 in each iteration. If the step is not accepted, then no update to the MINLP solution guess is performed, i.e., z^{k+1}=z^{k}, y^{k+1}=y^{k }366 and the trustregion radius value d_{k }can be updated 345 accordingly. Alternatively, if the step is accepted, then an update to the MINLP solution guess is performed for both integer and continuous variables, i.e., z^{k+1}=z^{k}+Δz^{k}, y^{k+1}=y+Δy^{k }365 and the trustregion radius value d_{k }can be updated 345 if necessary. In some embodiments of the invention, the check to accept a step or not 355 can be based on a merit function that quantifies a combination of optimality and constraint satisfaction for a solution guess of values for continuous and discrete decision variables. For example, if an update of the MINLP solution guess (Δy^{k},Δz^{k}) results in a sufficient decrease in the value of a merit function, then the step is accepted 355 and an update to the MINLP solution guess is performed 365, otherwise, if the decrease in the value of a merit function is insufficiently large, then the step is not accepted 355 and no update is performed 366.
Subsequently, the second iteration includes a partial convexification step within a subregion of the MINLP solution space that corresponds to a local neighborhood of (y^{1},z^{1}) 402, followed by a solution of the MICP subproblem to compute a search direction (Δy^{1},Δz^{1}) 412, a resulting update of the solution guess for integer variables z^{2}=z^{1}+Δz^{1 }422 and a resulting update of the solution guess for continuous variables y^{2}=y^{1}+α^{1}Δy^{1 }432. Similarly, one or multiple additional iterations can be performed until a termination condition 350 is satisfied, e.g., this condition may include checking whether a current intermediate solution guess is feasible, whether it is sufficiently close to the globally optimal solution and/or whether the computational cost of the sequential convexificationbased optimization procedure has reached a particular time limit.
Some embodiments of the invention are based on the realization that the solution of a convex QP subproblem 453 is generally much computationally cheaper than the solution of a nonconvex MIQP subproblem 451452 in each iteration of the sequential convexificationbased optimization procedure. In some embodiments of the invention, the decision whether to fix the values for all the integer variables 460 is based on whether the current intermediate solution guess is sufficiently close to the globally optimal solution and/or whether the computational cost of the iterative optimization procedure has reached a particular time limit. For example, in some embodiments, a maximum number N_{micp }of MICP subproblem solutions can be imposed to considerably reduce the computational effort of the iterative optimization procedure, e.g., by fixing the values of the integer variables Δz^{k}=0 after k≥N_{miqp }iterations 460.
For example, in
In some embodiments of the invention, each iteration of a sequential convexificationbased optimization procedure performs a partial convexification step to compute a convex approximation 525 of one or multiple nonconvex constraints 510 in a local neighborhood of a current solution guess 520, in order to compute a search direction that is approximately inside 526 the nonconvex feasible region 511. Some embodiments of the invention are based on the realization that each iteration of a sequential convexificationbased optimization procedure computes a solution to an MICP subproblem 320 that is nonconvex due to one or multiple integer variables, for example, including a binary decision variable z_{k}∈{0,1} 505, which causes a convex approximation 525 to exist for a value z_{k}=0 506, while a different convex approximation 535 of one or multiple nonconvex constraints 515 can exist in a local neighborhood of a current solution guess 530 corresponding to a value z_{k}=1 507, in order to compute a search direction that is approximately inside 536 the nonconvex feasible region 516.
For example, a solution to the MICP subproblem 320 can include setting a binary variable z_{k}∈{0,1} 505 to either z_{k}=0 506 or z_{k}=1 507, given the corresponding convex approximations 525 or 535 of a smooth nonlinear part of the MINLP. Some embodiments of the invention are based on the realization that a complex transformation 540 may exist between a current solution guess 520 or 530, a nonconvex feasible region 511 or 516, and a corresponding convex approximation 525 or 535 for each iteration of a sequential convexificationbased optimization procedure, depending on the value of a binary decision variable z_{k}∈{0,1} 505, i.e., depending on whether z_{k}=0 506 or z_{k}=1 507.
In some embodiments of the invention, a partial convexification step in each iteration of a sequential convexificationbased optimization procedure is based on a local linearization of one or multiple smooth nonlinear inequality constraint functions that define a nonconvex set of feasible values 511 or 516. Specifically, a partial convexification step can compute one or multiple linear inequality constraints, i.e., to stay inside one or multiple halfspaces 526 with respect to one or multiple linear functions 525, in order to approximate a nonconvex set of feasible values 511, defined by one or multiple smooth nonlinear inequality constraints 284285, within a subregion of the MINLP solution space that corresponds to a local neighborhood of the current intermediate solution guess 520. Similarly, a partial convexification step can compute one or multiple linear equality constraints in order to approximate a nonconvex set of feasible values, defined by one or multiple smooth nonlinear equality constraints 283, within a subregion of the MINLP solution space that corresponds to a local neighborhood of the current intermediate solution guess.
In some embodiments of the invention, each iteration of a sequential convexificationbased optimization procedure performs a partial convexification step to compute a convex approximation 565 of one or multiple nonconvex cost functions 550 in a local neighborhood of a current solution guess 560, in order to compute a search direction that is approximately towards a local 561 or global minimum 562 of the nonconvex cost function 550. Some embodiments of the invention are based on the realization that each iteration of a sequential convexificationbased optimization procedure computes a solution to an MICP subproblem 320 that is nonconvex due to one or multiple integer variables, for example, including a binary decision variable z_{k}∈{0,1} 505, which causes a convex approximation 565 to exist for a value z_{k}=0 506, while a different convex approximation 575 of one or multiple nonconvex cost functions 555 can exist in a local neighborhood of a current solution guess 570 corresponding to a value z_{k}=1 507, in order to compute a search direction that is approximately towards a local 571 or global minimum 572 of the nonconvex cost function 545.
Some embodiments of the invention are based on the realization that a complex transformation 580 may exist between a current solution guess 560 or 570, a nonconvex cost function 550 or 555, and a corresponding convex approximation 565 or 575 for each iteration of a sequential convexificationbased optimization procedure, depending on the value of a binary decision variable z_{k}∈{0,1} 505, i.e., depending on whether z_{k}=0 506 or z_{k}=1 507.
In some embodiments of the invention, a partial convexification step in each iteration of a sequential convexificationbased optimization procedure is based on a local linear or linearquadratic approximation 565 or 575 of one or multiple smooth nonlinear objective functions 281 that define a nonconvex cost function 550 or 555. Specifically, a partial convexification step can compute one or multiple linear or linearquadratic objective functions 565 or 575 in order to approximate a nonconvex cost function 545 within a subregion of the MINLP solution space that corresponds to a local neighborhood of the current intermediate solution guess 560 or 570, and the resulting MICP subproblem includes the one or multiple linear or linearquadratic objective functions 565 or 575 such that the solution of the MICP subproblem defines a search direction within the MINLP solution space towards a local or global minimum of a nonconvex cost function 545.
Based on a compact formulation y=[x_{0}^{T}, u_{0}^{T}, x_{1}^{T}, . . . , u_{N−1}^{T}, x_{N}^{T}]^{T }of continuous variables 636 and z=[w_{0}^{T}, w_{1}^{T}, . . . , w_{N}^{T}]^{T }of integer variables 637, a compact formulation of the optimal control structured MINLP in separable format 630 can be formulated, and which is mathematically equivalent 625 to the original optimal control structured MINLP in separable format 280, i.e., the compact MINLP 630 is infeasible if and only if the original MINLP 280 is infeasible, and a feasible and (locally) optimal solution to the compact MINLP 630 can be used to construct a feasible and (locally) optimal solution of the original MINLP 280.
In some embodiments of the invention, the objective function 631 in the compact formulation of the MINLP 630 is defined by a smooth nonlinear function ƒ(y): ^{n}^{y}→, which depends only on the continuous variables y∈^{n}^{y }636, and a linear function c^{T }z that is defined b a vector c∈^{n}^{z }corresponding to the integer variables z∈^{n}^{z }637. Similarly, the equality constraints 632 in the compact formulation of the MINLP 630 can be defined by a smooth nonlinear function g(y): ^{n}^{y}→^{n}^{g }and a set of linear equations Dz that is defined by a matrix D∈^{n}^{g}^{×n}^{z }corresponding to the integer variables z∈^{n}^{z }637. Similarly, the inequality constraints 633 in the compact formulation of the MINLP 630 can be defined by a smooth nonlinear function h(y): ^{n}^{y}→^{n}^{h }and a set of linear equations Ez that is defined by a matrix E∈^{n}^{h}^{×n}^{z }corresponding to the integer variables z∈^{n}^{z }637. Finally, the compact formulation of the MINLP 630 can include one or multiple discrete equality constraints 634, which impose that each variable z_{j}∈ can only take integer values in the feasible MINLP solution space.
In some embodiments of the invention, the MIQP subproblem 650 is constructed based on a sequential convex programming (SCP) method or a sequential quadratic programming (SQP) method applied to the smooth nonlinear functions in the compact formulation of the MINLP in separable format 630 with respect to the continuous variables y 636, i.e., without affecting the linear terms with respect to the integer variables z 637 of the MINLP. Some embodiments of the invention are based on the realization that the MINLP in separable format 630 corresponds to a smooth nonlinear programming (NLP) problem, when fixing the integer variables to a feasible and/or (locally) optimal set of values z=
In some embodiments of the invention, in each iteration of the sequential convexificationbased optimization procedure, the objective 651 of the MIQP subproblem 650 consists of a linearquadratic term that depends on a gradient vector ∇_{y}ƒ(y^{k}) for the function ƒ(y): ^{n}^{y}→, evaluated at the current solution guess y^{k}∈^{n}^{y}, and on a symmetric Hessian matrix B(y^{k}) that is positive semidefinite in general, or at least positive definite in the null space of all equality and all active inequality constraints to ensure that the MIQP subproblem 650 becomes a convex QP for a fixed set of values for integer variables Δz=
(y,λ,μ)=ƒ(y)+λ^{T}g(y)+μ^{T}h(y),
where λ denote the Lagrange multipliers with respect to the equality constraints 632 and μ denote the Lagrange multipliers with respect to the inequality constraints 633 of the MINLP in separable format 630.
In some embodiments of the invention, a possible regularization term can be added to the Hessian matrix B(y^{k}) to ensure that the resulting matrix is positive semidefinite, for example, B(y^{k})+γ0, where γ≥0 is a nonnegative regularization value and is an identity matrix. In some embodiments of the invention, a quasiNewton Hessian approximation method can be used to compute a computationally efficient approximation to the Hessian of the Lagrangian, e.g., based on a symmetric rank1 (SR1) update method or using a variant of the BroydenFletcherGoldfarbShanno (BFGS) method. In some embodiments of the invention, e.g., if the objective 631 consists of a nonlinear least squares function ƒ(y)=∥r(y)∥_{2}^{2}, then a GaussNewton Hessian approximation to the Hessian of the Lagrangian can be computed as
denotes the Jacobian of a function r(y): ^{n}^{y}→^{n}^{r }with respect to the continuous variables y 636 and evaluated at a current solution guess y^{k}∈^{n}^{y}. Some embodiments of the invention are based on the realization that a GaussNewton Hessian approximation is computationally cheap to evaluate and is guaranteed to be positive semidefinite, i.e., B(y^{k})0.
In some embodiments of the invention, the MIQP subproblem 650 at each iteration of the optimization procedure includes one or multiple linear equality constraints 652, based on a local linearization of the nonlinear equality constraints 632 of the MINLP in separable format 630. Specifically, the smooth nonlinear function g(y): ^{n}^{y}→^{n}^{g }can be approximated by a linearization
in a local neighborhood around a current solution guess y^{k}∈^{n}^{y}, where
denotes the Jacobian of a function g(y): ^{n}^{y}→^{n}^{g }with respect to the continuous variables y 636 and evaluated at a current solution guess y^{k}∈^{n}^{y}. In addition, the linear term Dz in the nonlinear equality constraints 632 remains unchanged in the linearized equality constraints 652, i.e., Dz^{k+1}=D(z^{k}+Δz).
In some embodiments of the invention, the MIQP subproblem 650 at each iteration of the optimization procedure includes one or multiple linear inequality constraints 653, based on a local linearization of the nonlinear inequality constraints 633 of the MINLP in separable format 630. Specifically, the smooth nonlinear function h(y): ^{n}^{y}→^{n}^{h }can be approximated by a linearization
in a local neighborhood around a current solution guess y^{k}∈^{n}^{y}, where
denotes the Jacobian of a function h(y): ^{n}^{y}→^{n}^{h }with respect to the continuous variables y 636 and evaluated at a current solution guess y^{k}∈^{n}^{y}. In addition, the linear term Ez in the nonlinear inequality constraints 633 remains unchanged in the linearized inequality constraints 653, i.e., Ez^{k+1}=E(z^{k}+Δz). Similarly, the discrete equality constraints 634 of the MINLP in separable format 630 remain unchanged in the discrete equality constraints 654 of the MIQP subproblem 650.
Some embodiments of the invention are based on the realization that the solution of a convex QP subproblem 660 is generally much computationally cheaper than the solution of a nonconvex MIQP subproblem 650 in each iteration of the sequential convexificationbased optimization procedure. In some embodiments of the invention, the decision whether to fix the values for all the integer variables 665 is based on whether the current intermediate solution guess y^{k}∈^{n}^{y}, z^{k}∈^{n}^{z }is sufficiently close to the globally optimal solution and/or whether the computational cost of the iterative optimization procedure has reached a particular time limit. For example, in some embodiments, a maximum number N_{miqp }of MIQP subproblem solutions can be imposed to considerably reduce the computational effort of the iterative optimization procedure, e.g., by fixing the values of the integer variables Δz^{k}=0 after k≥N_{miqp }iterations 665.
In some embodiments of the invention, a sequential convexificationbased optimization procedure solves a convex QP subproblem 660 in one or multiple subsequent iterations to update the MINLP solution guess y^{k+1}∈^{n}^{y}, z^{k+1}=z^{k}∈^{n}^{z }641 after fixing the values for integer variables 665, until a feasible and (locally) optimal solution has been found 255 for the MINLP in separable format. In some embodiments of the invention, the integer variables can be turned back into free optimization variables, after one or multiple iterations of a sequential convexificationbased optimization procedure based on a convex QP subproblem solution 660, for example, if the MINLP is detected to be infeasible for a fixed set of values z=
If the current solution guess for the integer variables z^{k}∈^{n}^{z }is sufficiently good 670, then the values for integer variables can be fixed, i.e., Δz^{k}=0 and a convex QP subproblem can be constructed and solved to compute Δy^{k }675. Alternatively, if the current solution guess for the integer variables z^{k}∈^{n}^{z }is not sufficiently good 670, then a nonconvex MIQP subproblem is constructed and solved to compute (Δy^{k},Δz^{k}) 680. Given a search direction (Δy^{k},Δz^{k}), computed by either solving a convex QP subproblem 675 or by solving a nonconvex MIQP subproblem 680, it can be used directly to update the intermediate solution guess for integer variables z^{k+1}=z^{k}+Δz^{k }325, and it can be used to select a step size 0≤α^{k}≤1 to update the intermediate solution guess for continuous variables y^{k+1}=y^{k}+α^{k}Δy^{k }330, resulting in an updated MINLP solution guess y^{k+1}∈^{n}^{y}, z^{k+1}∈^{n}^{z }335 that can be used in the next iteration k←k+1 340.
F(y,z)=ƒ(y)+c^{T}z,
G(y,z)=g(y)+Dz, H(y,z)=h(y)+Ez,
a merit function can be defined as follows 711
where a value of the merit function can be minimized to find a feasible and (locally) optimal solution of the MINLP in separable format 630.
In some embodiments of the invention, the merit function can consist of one or multiple terms. For example, a first term in the merit function 711 can correspond to the objective function F(y,z)=ƒ(y)+c^{T}z 631, which quantifies the optimality of an MINLP solution guess (y,z) in the MISCP optimization algorithm. A second term in the merit function 711 can be defined as ρ∥G(y,z)∥_{1}, where ρ>0 denotes a positive penalty parameter value and G(y,z)=0 is a compact notation 705 for the nonlinear equality constraints g(y)+Dz=0 632, such that a term ρ∥G(y,z)∥_{1 }quantifies the satisfaction of the nonlinear equality constraints for an MINLP solution guess (y,z) in the MISCP optimization algorithm. If an MINLP solution guess is feasible, then G(y,z)=g(y)+Dz=0 and therefore the second term is minimized and ρ∥G(y,z)∥_{1}=0. In some embodiments of the invention, a different norm can be used in the merit function such as, for example, a 2norm or an ∞norm instead of the 1norm.
Finally, a third term in the merit function 711 can be defined as ρΣ_{i }max(H_{i}(y,z),ϵ), where ϵ≥0 is a small nonnegative tolerance value and H(y,z)≤0 is a compact notation 705 for the nonlinear inequality constraints h(y)+Ez≤0 633, such that a term ρΣ_{i }max(H_{i}(y,z),ϵ) quantifies the satisfaction of the nonlinear inequality constraints for an MINLP solution guess (y,z) in the MISCP optimization algorithm. If an MINLP solution guess is feasible, then H_{i}(y,z)=h_{i}(y)+E_{i}, z≤0 and therefore the third term is minimized and ρΣ_{i }max(H_{i}(y,z),ϵ)=ρΣ_{i }ϵ. In some embodiments of the invention, the tolerance value ϵ=0 but, in other embodiments, the tolerance value can be chosen as a small positive value, for example, ϵ=10^{−6}.
In some embodiments of the invention, the MISCP optimization algorithm performs one or multiple iterations, in which each iteration performs a partial convexification step, and then constructs and solves an MICP subproblem 715 to compute a search direction for continuous variables Δy^{k }and a search direction for integer variables Δz^{k}. In some embodiments of the invention, given the solution of the MICP subproblem 715, the solution guess for integer variables can be updated directly as z^{k+1}=z^{k}+Δz^{k}. Some embodiments of the invention are based on the realization that a search direction for continuous variables Δy^{k}, computed as the solution of an MICP subproblem 715 such as the MIQP subproblem 650 in
∇_{y}ϕ(y^{k};z^{k+1},ρ)^{T}Δy^{k}<0.
given a sufficiently large parameter value ρ>0 720.
In some embodiments of the invention, each iteration of the MISCP optimization algorithm computes a sufficiently large parameter value ρ>0 such that the descent condition 721 holds. Embodiments of the invention are based on the realization that increasing the parameter value ρ>0 results in a merit function 711 that quantifies constraint satisfaction relatively more strongly compared to optimality for the MINLP solution guess. Alternatively, decreasing the parameter value ρ>0 results in a merit function 711 that quantifies optimality relatively more strongly compared to constraint satisfaction for the MINLP solution guess.
In some embodiments of the invention, the step size value 0≤α^{k}≤1 is selected to ensure a sufficient decrease condition holds for a merit function of the MINLP, for example the 1norm based merit function 711, as follows 731
ϕ(y^{k}+α^{k}Δy^{k},z^{k+1},ρ)≤ϕ(y^{k},z^{k+1},ρ)+α^{k}η∇_{y}ϕ(y^{k};z^{k+1},ρ)^{T}Δy^{k},η∈(0,1),
where a parameter value 0<η<1 can be chosen. Embodiments of the invention are based on the realization that the sufficient decrease condition 731 ensures that a merit function value ϕ(y^{k}+α^{k}Δy^{k};z^{k+1},ρ) for the updated MINLP solution guess y^{k+1}=y^{k}+α^{k}Δy^{k }730 is at least smaller than the merit function value ϕ(y^{k};z^{k+1},ρ) for the current MINLP solution guess y^{k}, given the directional derivative computation ∇_{y}ϕ(y^{k};z^{k+1},ρ) for the merit function evaluated at the current solution guess y^{k }for continuous variables and the updated solution guess z^{k+1 }for integer variables. Embodiments of the invention are based on the realization that the directional derivative of a merit function can be evaluated computationally efficiently using symbolic differentiation, algorithmic differentiation, and/or numerical differentiation tools. For example, the gradient ∇_{y}ϕ(y^{k};z^{k+1},ρ) can be evaluated efficiently using an adjoint mode of algorithmic differentiation applied to the merit function 711.
In some embodiments of the invention, an iterative backtracking procedure is used to select the step size value 0≤α^{k}≤1, starting from an initial value α^{k}=1 and this step size value can be decreased iteratively towards 0≤α^{k }until the sufficient decrease condition 731 holds for a merit function of the MINLP.
In some embodiments of the invention, the termination condition of the MISQP optimization algorithm can be based on a norm of the KarushKuhnTucker (KKT) necessary conditions of optimality for the MINLP ∥r(y^{k},z^{k})∥ 742, excluding the integrality conditions. Therefore, as long as the condition ∥r(y^{k},z^{k})∥>ϵ_{tol }742 holds, the MISQP optimization algorithm performs one or multiple additional iterations to compute a feasible and (locally) optimal MINLP solution.
Each iteration performs a partial convexification step 315, followed by the solution of an MIQP subproblem to compute a search direction (Δy^{k},Δz^{k}) 680 such as the MIQP subproblem 650 in
Given a new search direction (Δy^{k},Δz^{k}), each iteration of the MISQP optimization method can update the solution guess for integer variables z^{k+1}=z^{k}+Δz^{k }325, an iteration can compute a sufficiently large parameter value ρ>0 such that a descent condition 721 holds, and a step size value can be initialized as α^{k}←1 747. Then, in some embodiments of the invention, an iterative backtracking procedure is used to select the step size value 0≤α^{k}≤1 while a sufficient decrease condition 731 is not yet satisfied 748. For example, each iteration of the iterative backtracking procedure decreases the current step size value α^{k}←{tilde over (β)}α^{k }749, given a parameter value 0<{tilde over (β)}≤β<1, until a sufficient decrease condition 731 is satisfied for a merit function of the MINLP. In some embodiments of the invention, the value for {tilde over (β)} is decreased gradually in each iteration of the iterative backtracking procedure in order to reduce the number of evaluations for the merit function of the MINLP, and therefore reduce the computational cost of the step size selection in the MISQP optimization algorithm. At the end of each iteration, the step size value 0≤α^{k}≤1 can be used to update the solution guess for continuous variables y^{k+1}=y^{k}+α^{k}Δy^{k }730 to obtain a new MINLP solution guess y^{k+1}∈^{n}^{y}, z^{k+1}∈^{n}^{z }335 for the next iteration of the MISQP optimization algorithm 740.
where a value of the merit function can be minimized to find a feasible and (locally) optimal solution of the MINLP in separable format 630.
In some embodiments of the invention, a linearizationbased approximation of the merit function 805 can be used, given the MINLP solution guess (y^{k},z^{k}) and MICP subproblem solution (Δy^{k},Δz^{k}), to predict the optimality and constraint satisfaction in the search direction that is computed by the MICP subproblem solution, for example, as follows
which includes a linearquadratic approximation of the objective function 806, based on a gradient vector ∇_{y}ƒ(y^{k}) for the function ƒ(y): ^{n}^{y}→, evaluated at the current solution guess y^{k}∈^{n}^{y}, and on a symmetric Hessian matrix B(y^{k}) that is positive semidefinite in general, similar to the linearquadratic objective 651 of the MIQP subproblem 650 in
652 of the MIQP subproblem 650 in
653 of the MIQP subproblem 650 in
In some embodiments of the invention, given a merit function for the MINLP 711, given a linearizationbased approximation of the merit function 805, given an MINLP solution guess y^{k}∈^{n}^{y}, z^{k}∈^{n}^{z }335, given a search direction Δy^{k}∈^{n}^{y }and given an updated solution guess for integer variables z^{k+1}=z^{k}+Δz^{k }815, a ratio R_{k }of actual to predicted reduction for the value of a merit function 810 can be defined, for example, as follows 811
where a positive ratio value R_{k}>0 can correspond to a reduction in the actual value of the merit function ϕ(y^{k}+Δy^{k};z^{k+1},ρ)<ϕ(y^{k};z^{k+1},ρ) 710 and a reduction in the linearizationbased approximation of the merit function ϕ_{QP}^{k}(Δy^{k};z^{k+1},ρ)<ϕ_{QP}^{k}(0;z^{k+1},ρ) 805. Some embodiments of the invention are based on the realization that a value for the ratio R_{k }needs to be positive and sufficiently large in order to accept a step in the search direction Δy^{k }for the continuous variables in each iteration of the MISCP optimization method.
In some embodiments of the invention, the update of a new trustregion radius value d_{k+1 }820 can include decreasing the trustregion radius value d_{k+1 }820, compared to the current radius value d_{k}, if the ratio value R_{k}<η_{1}<<1 861 is considerably small or even negative. In addition, in some embodiments of the invention, the update of a new trustregion radius value d_{k+1 }820 can include decreasing the trustregion radius value d_{k+1 }820, compared to the current radius value d_{k}, if the step search direction Δy^{k }for the continuous variables is considerably smaller than the current radius value d_{k}, for example, if ∥MΔy_{k}∥_{p}<γ_{1}d_{k }863 where M0 denotes a scaling matrix and 0<γ_{1}<1. In some embodiments of the invention, the value p∈[1,∞] defines the norm for the trustregion procedure in the MISCP optimization algorithm, for example, p=1 or p=∞ to define the 1norm or ∞norm, respectively.
In some embodiments of the invention, the update of a new trustregion radius value d_{k+1 }820 can include increasing the trustregion radius value d_{k+1 }820, compared to the current radius value d_{k}, for example, if the ratio value R_{k }is sufficiently large and positive R_{k}>η_{2}>0 and the step search direction Δy^{k }for the continuous variables is strictly bounded by the current trustregion radius value d_{k}, i.e., the condition ∥MΔy^{k}∥_{p}=d_{k }865 holds as an equality constraint.
Some embodiments of the invention are based on the realization that, in general, the MICP subproblem becomes computationally cheaper to solve for increasingly small values of the trustregion radius d_{k}. Therefore, in some embodiments of the invention, the trustregion procedure aims to reduce the trustregion radius when it is possible without slowing down convergence of the MISCP optimization algorithm to a feasible and (locally) optimal MINLP solution.
Each iteration of a trustregion MISQP optimization method performs a partial convexification step 315, followed by the solution of an MIQP subproblem to compute a search direction (Δy^{k},Δz^{k}) 855 such as the MIQP subproblem 840 in
Given a new search direction (Δy^{k},Δz^{k}), each iteration of the trustregion MISQP optimization method can compute a ratio R_{k }of actual to predicted reduction for the value of a merit function 810 and, if the ratio value R_{k }is sufficiently large and positive R_{k}≥η_{1}>0 825 then the search direction is accepted and a fullstep update to the MINLP solution guess y^{k+1}←y^{k}+Δy^{k }and z^{k+1}←z^{k}+Δz^{k }830 is computed and otherwise, if the ratio value R_{k }is not sufficiently large and positive R_{k}<η_{1 }825 then the search direction is rejected and the MINLP solution guess is not updated, i.e., y^{k+1}←y^{k }and z^{k+1}←z^{k }831. Next, in some embodiments of the invention, a new trustregion radius value is computed 820.
Some embodiments of the invention shrink the trustregion radius value d_{k+1}←max(γ_{1}∥MΔy^{kl}∥_{p},d) 862, if the ratio value R_{k}<η_{1}<<1 861 is considerably small or even negative. In addition, some embodiments of the invention shrink the trustregion radius value d_{k+1}←max(γ_{2}d_{k},d) 864, if the step search direction Δy^{k }for the continuous variables is considerably smaller than the current radius value d_{k}, for example, if ∥MΔy^{k}∥_{p}<γ_{1}d_{k }863 where M0 denotes a scaling matrix and 0<γ_{1}<1. Some embodiments of the invention, for example, grow the trustregion radius value d_{k+1}←min(γ_{3}d_{k},
Some embodiments of the invention are based on the realization that enforcing an optimization variable to be binary, i.e., z_{j}∈{0,1} 900, can be equivalent 905 to adding a penalty term β^{k}z_{j}(1−z_{j}) 911 to the MINLP objective function 631, while restricting 0≤z_{j}≤1 912, for a homotopy sequence of increasing parameter values β^{k+1}≥β^{k }such that β^{k}→∞ for iterations k→∞ 910. For example, in some embodiments of the invention, the MISCP optimization algorithm performs one or multiple iterations including a penalty term β^{0}z_{j}(1−z_{j}) 914 in the MINLP objective function 631, followed by an update to the homotopy penalty parameter value β^{1}≥β^{0 }915, then the MISCP optimization algorithm performs one or multiple iterations including a penalty term β^{1}z_{j}(1−z_{j}) 916 in the MINLP objective function 631, followed by an update to the homotopy penalty parameter value β^{2}≥β^{1 }917, and these computational steps can be repeated 918919 for one or multiple iterations of the MISCP optimization algorithm until a feasible and (locally) optimal solution is found for the MINLP.
In addition, in some embodiments of the invention, one or multiple linear terms w_{j}(1−z_{j}) 940 can be added to the cost function of the MICP minimization subproblem, where w_{j}>0 is a positive weight value, to ensure a binary variable z_{j }921, which is bounded as 0≤z_{j}≤1 912, to remain close to a solution guess z_{j}^{k}=1. In some embodiments of the invention, a linear term w_{j}^{0}(1−z_{j}) 940 in the cost function of the MICP minimization subproblem can be computed as a linear approximation of a penalty term β^{0}z_{j}(1−z_{j}) 914 at an MINLP solution guess 932. Similarly, after one or multiple iterations of the MISCP optimization algorithm, and for an increased homotopy penalty parameter value β^{1}≥β^{0 }915, a linear term w_{j}^{1}(1−z_{j}) 941 in the cost function of the MICP minimization subproblem can be computed as a linear approximation of a penalty term β^{1}z_{j}(1−z_{j}) 916 at an updated MINLP solution guess 933.
Some embodiments of the invention are based on the realization that the use of a homotopytype penalty method for some or all of the integer and/or binary decision variables may prevent cycling in the MISCP optimization algorithm. In addition, some embodiments of the invention are based on the realization that the use of a homotopytype penalty method for some or all of the integer and/or binary decision variables may considerably reduce the computational cost of solving the MICP subproblems in the MISCP optimization algorithm. For example, some embodiments of the invention are based on a branchandbound method to solve the MICP subproblem in each iteration of the MISCP optimization algorithm, and adding one or multiple linear penalty terms to the cost function may result in a considerably smaller branchandbound search tree and therefore a considerably reduced computational cost of solving the MICP subproblem.
Some embodiments of the invention are based on the realization that enforcing an optimization variable to be binary, i.e., z_{j}∈{0,1} 900, can be equivalent 955 to adding a smooth nonlinear inequality constraint β^{k}z_{j}(1−z_{j})≤1 951 to the MINLP constraints 633, while restricting 0≤z_{j}≤1 912, for a homotopy sequence of increasing parameter values β^{k+1}≥β^{k }such that β^{k}→∞ for iterations k→∞ 950. For example, in some embodiments of the invention, the MISCP optimization algorithm performs one or multiple iterations including a smooth nonlinear inequality constraint β^{0}z_{j}(1−z_{j})≤1 961 in the MINLP constraints 633, followed by an update to the homotopy penalty parameter value β^{1}≥β^{0 }962, then the MISCP optimization algorithm performs one or multiple iterations including a smooth nonlinear inequality constraint β^{1}z_{j}(1−z_{j})≤1 963 in the MINLP constraints 633, followed by an update to the homotopy penalty parameter value β^{2}≥β^{1 }964, and these computational steps can be repeated 965966 for one or multiple iterations of the MISCP optimization algorithm until a feasible and (locally) optimal solution is found for the MINLP.
In some embodiments of the invention, a homotopytype penalty method can be used to add one or multiple inequality constraints to each MICP subproblem 320 to increasingly enforce the MISCP optimization algorithm to compute an update to the MINLP solution guess for some or all of the integer and/or binary decision variables that remains close to the solution guess of the integer and/or binary decision variables at a previous MISCP iteration. For example, a smooth nonlinear inequality constraint β^{0}z_{j}(1−z_{j})≤1 970 may be feasible for all values of 0≤z_{j}≤1 912, given a particular homotopy parameter value β^{0}>0. However, for an increasingly large homotopy penalty parameter value β^{1}≥β^{0 }962, a subregion of the solution space 975 may become infeasible, and an increasingly large subregion of the solution space 976 may become infeasible for an increasingly large homotopy penalty parameter value β^{2}≥β^{1 }964, such that the MINLP solution guess is enforced to become equal to either z_{j}^{k}=0 or z_{j}^{k}=1 in one or multiple iterations of the MISCP optimization algorithm.
In some embodiments of the invention, if a particular set of values for one or multiple integer variables and/or continuous variables is detected to be infeasible for the MINLP 630, then each of the subsequent iterations of the MISCP optimization algorithm can include one or multiple linear and/or smooth nonlinear inequality constraints to avoid revisiting the same set of values for one or multiple integer variables and/or continuous variables for the MINLP 630.
Some embodiments of the invention are based on a branchandbound method to solve the MICP subproblem in each iteration of the MISCP optimization algorithm, and adding one or multiple inequality constraints may result in a considerably smaller branchandbound search tree and therefore a considerably reduced computational cost of solving the MICP subproblem.
{tilde over (y)}=[x_{0}^{T},u_{0}^{T},x_{1}^{T},u_{1}^{T}, . . . ,u_{N−1}^{T},x_{N}^{T}]^{T}, {tilde over (z)}=[w_{0}^{T},w_{1}^{T}, . . . ,w_{N−1}^{T}w_{N}^{T}]^{T }
which can be shifted in time by one sampling time period 1015 to result in the following shifted MINLP solution guess 1017
ŷ=[x_{1}^{T},u_{1}^{T}, . . . ,u_{N−1}^{T},x_{N}^{T},ϕ_{u}^{+}(u_{N−1})^{T},ϕ_{x}^{+}(x_{N},u_{N−1})^{T}]^{T},
{circumflex over (z)}=[w_{1}^{T}, . . . ,w_{N−1}^{T},w_{N}^{T},ϕ_{w}^{+}(w_{N})^{T}]^{T }
where ϕ_{u}^{+}(⋅), ϕ_{x}^{+}(⋅) and ϕ_{w}^{+}(⋅) denote linear or nonlinear operators to predict values for the control input variables, the state variables and the integer decision variables at the next time step within the prediction time horizon. In some embodiments of the invention, the time shifted set of optimal control values can be used directly as an initial MINLP solution guess (y^{0},z^{0})=(ŷ,{circumflex over (z)}) 1020 for the MISCP optimization algorithm to compute a feasible and (locally) optimal solution for the MINMPC problem at the current control time step t_{1}.
For example, the partition P_{1 }1101 represents a discrete search region that can be split or branched into two smaller partitions or regions P_{2 }1102 and P_{3 }1103, i.e., a first and a second region that are nested in a common region. The first and the second region are disjoint, i.e., the intersection of these regions is empty P_{2}∩P_{3}=ϕ 1107, but they form the original partition or region P_{1 }together, i.e., the union P_{2}∪P_{3}=P_{1 }1106 holds after branching. The branchandbound method then solves an integerrelaxed convex problem for both the first and the second partition or region of the search space, resulting in two solutions (local optimal solutions) that can be compared against each other as well as against the currently known upper bound value to the optimal objective value. The first and/or the second partition or region can be pruned if their performance metric is less optimal than the currently known upper bound to the optimal objective value of the MICP subproblem. The upper bound value can be updated if the first region, the second region or both regions result in a discrete feasible solution to the MICP subproblem. The branchandbound method then continues by selecting a remaining region in the current nested tree of regions for further partitioning.
While solving each partition may still be challenging, it is fairly efficient to obtain local lower bounds on the optimal objective value, by solving local relaxations of the mixedinteger program or by using duality. If the solution method for the MICP subproblem happens to obtain an integerfeasible solution while solving a local relaxation, the branchandbound method can then use it to obtain a global upper bound for the mixedinteger control solution of the original MICP subproblem in the MISCP optimization algorithm to find a feasible and (locally) optimal solution of the MINLP. This may help to avoid solving or branching certain partitions that were already created, i.e., these partitions or nodes can be pruned. This general algorithmic idea of partitioning can be represented as a binary search tree 1100, including a root node, e.g., P_{1 }1101 at the top of the tree, and leaf nodes, e.g., P_{4 }1104 and P_{5 }1105 at the bottom of the tree. In addition, the nodes P_{2 }1102 and P_{3 }1103 are typically referred to as the direct children of node P_{1 }1101, while node P_{1 }1101 is referred to as the parent of nodes P_{2 }1102 and P_{3 }1103. Similarly, nodes P_{4 }1104 and P_{5 }1105 are children of their parent node P_{2 }1102.
As long as the gap between the lower and upper bound value is larger than a particular tolerance value at step 1111, and a maximum execution time is not yet reached by the optimization algorithm, then the branchandbound method continues to search iteratively for the mixedinteger optimal control solution of the MICP subproblem 1155. Each iteration of the branchandbound method starts by selecting the next node in the tree, corresponding to the next region or partition of the integer variable search space, with possible variable fixings based on presolve branching techniques 1115. After the node selection, the corresponding integerrelaxed convex problem is solved, with possible variable fixings based on postsolve branching techniques 1120.
If the integerrelaxed convex problem has a feasible solution, then the resulting relaxed control solution provides a lower bound on the objective value for that particular region or partition of the integer variable search space. At step 1121, if the objective is determined to be larger than the currently known upper bound for the objective value of the optimal mixedinteger control solution of the MICP subproblem, then the selected node is pruned or removed from the branching tree 1140. However, at step 1121, if the objective is determined to be lower than the currently known upper bound, and the relaxed control solution is integer feasible 1125, then the currently known upper bound and corresponding mixedinteger control solution estimate is updated at step 1130.
If the integerrelaxed convex problem has a feasible solution and the objective is lower than the currently known upper bound 1121, but the relaxed control solution is not yet integer feasible 1125, then the global lower bound for the objective can be updated 1135 to be the minimum of the objective values for the remaining leaf nodes in the branching tree and the selected node is pruned from the tree 1140. In addition, starting from the current node, a discrete variable with a fractional value is selected for branching according to a particular branching strategy 1145, in order to create and append the resulting auxiliary MICP subproblems, corresponding to subregions or partitions of the discrete search space of the original MICP subproblem, as children of that node in the branching tree 1150.
An important step in the branchandbound method is how to create the partitions, i.e., which node to select 1115 and which discrete variable to select for branching 1145. Some embodiments are based on branching one of the binary control variables with fractional values in the integerrelaxed convex problem solution. For example, if a particular binary control variable u_{i,k}∈{0,1} has a fractional value as part of the integerrelaxed convex problem solution, then some embodiments create two partitions of the MICP subproblem by adding, respectively, the equality constraint u_{i,k}=0 to one subproblem and the equality constraint u_{i,k}=1 to the other subproblem. Some embodiments are based on a reliability branching strategy for variable selection 1145, which aims to predict the future branching behavior based on information from previous branching decisions.
Some embodiments are based on a branchandbound method that uses a depthfirst node selection strategy, which can be implemented using a lastinfirstout (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, i.e., the node is either infeasible, optimal or dominated by the currently known upper bound value, which is followed by a backtracking procedure. Instead, some embodiments are based on a branchandbound method that uses a bestfirst strategy that selects the node with the currently lowest local lower bound. Some embodiments employ a combination of the depthfirst and bestfirst node selection approach, in which the depthfirst node selection strategy is used until an integerfeasible control solution is found, followed by using the bestfirst node selection strategy in the subsequent iterations of the branchandbound based optimization algorithm. The latter implementation is motivated by aiming to find an integerfeasible control solution early at the start of the branchandbound procedure (depthfirst) to allow for early pruning, followed by a more greedy search for better feasible solutions (bestfirst).
The branchandbound method continues iterating until either one or multiple of the following conditions have been satisfied:
The maximum execution time for the processor is reached.
All the nodes in the branching search tree have been pruned, such that no new node can be selected for solving convex relaxations or branching.
The optimality gap between the global lower and upper bound value for the objective of the MICP subproblem solution is smaller than the tolerance.
The vehicle can also include an engine 1206, which can be controlled by the controller 1202 or by other components of the vehicle 1201. The vehicle can also include one or more sensors 1204 to sense the surrounding environment. Examples of the sensors 1204 include distance range finders, radars, lidars, and cameras. The vehicle 1201 can also include one or more sensors 1205 to sense its current motion quantities and internal status. Examples of the sensors 1205 include global positioning system (GPS), accelerometers, inertial measurement units, gyroscopes, shaft rotational sensors, torque sensors, deflection sensors, pressure sensor, and flow sensors. The sensors provide information to the controller 1202. The vehicle can be equipped with a transceiver 1206 enabling communication capabilities of the controller 1202 through wired or wireless communication channels.
The spacecraft 1302 flies in outer space along an open or closed orbital path 1360 around, between, or near one or more gravitational bodies such as the Earth 1361, moon, and/or other celestial planets, stars, asteroids, comets. Usually, a desired or target position 1365 along the orbital path is given. A reference frame 1370 is attached to the desired position, where the origin of the frame, i.e., the all zeros coordinates in that reference frame are the coordinates of the desired position at all times.
The spacecraft 1302 is subject to various disturbance forces 1314. These disturbance forces can include forces that were not accounted for when determining the orbital path for the spacecraft 1302. These disturbance forces act on the spacecraft 1302 to move the spacecraft 1302 away from the desired position on the orbital path. These forces can include, but are not limited to, gravitational attraction, radiation pressure, atmospheric drag, nonspherical central bodies, and leaking propellant. Thus, the spacecraft 1302 can be at a distance 1367 away from the target position.
Because of the disturbance forces, it is not always possible to keep the spacecraft 1302 at the desired position along its orbit. As such, it is desired that the spacecraft 1302 instead remains within a window 1366 with specified dimensions 1364 around the desired position. To that end, the spacecraft 1302 is controlled to move along any path 1380 that is contained within the desired target window. In this example, the window 1366 has a rectangular shape, but the shape of the window can vary for different embodiments.
The spacecraft 1302 is also often required to maintain a desired orientation. For example, a spacecraftfixed reference frame 1374 is required to be aligned with a desired reference frame such as an inertial reference frame 1371 that is fixed relative to distant stars 1372, or a reference frame 1373 that is always oriented in a manner that points towards the Earth. However, depending on the shape of the spacecraft 1302, different disturbance forces 1314 can act nonuniformly on the spacecraft 1302, thereby generating disturbance torques, which cause the spacecraft 1302 to rotate away from its desired orientation. In order to compensate for the disturbance torques, momentum exchange devices 1351 such as reaction wheels are used to absorb the disturbance torques, thus allowing the spacecraft to maintain its desired orientation.
So that the momentum exchange devices do not saturate, and thereby lose the ability to compensate for disturbance torques, their stored momentum must be unloaded, e.g., by reducing spin rates of the reaction wheels. Unloading the momentum exchange devices imparts an undesired torque on the spacecraft 1302. Such an undesired torque is also compensated for by the thrusters.
In some embodiments of the invention, using one or multiple nonlinear equations to describe the dynamics and/or constraints of the controlled system, the MINMPC controller determines an input to the spacecraft 1302 based on the mixedinteger optimal control solution, wherein the input to the spacecraft 1302 actuates one or a combination of thrusters and momentum exchange devices, and the discrete optimization variables are used to model one or a combination of discrete control decisions, switching in the system dynamics, integer values for the thruster commands, and obstacle avoidance constraints.
In some embodiments of the invention, the spacecraft 1302 can be modeled as a hybrid nonlinear system and the commands that are sent to the actuators are computed using a predictive controller, such as the mixedinteger nonlinear model predictive controller (MINMPC). For example, in some embodiments, the commands that are sent to the thrusters 1350 can only take a discrete set of values, and therefore resulting into a set of binary or integer control input variables for each stage within the mixedinteger control horizon.
In some embodiments of the invention, the predictive controller is designed such that the spacecraft 1302 remains outside of a particular zone 1385 with specified dimensions, close to the desired position along the orbit. The latter zone can be either fixed in time or it can be time varying, and is often referred to as an exclusion zone 1385, for which the corresponding logic inequality constraints can be modeled using an additional set of binary or integer control input variables for each stage within the mixedinteger control horizon. In this example, the exclusion zone 1385 has a rectangular shape, and it is positioned in a corner of the desired window 1366, but the shape and position of the exclusion zone within the desired target window can vary for different embodiments.
Additionally, the VCS 1400 can include a flow reversing valve 1455 that is used to direct high pressure refrigerant exiting the compressor to either the outdoor unit heat exchanger or the indoor unit heat exchanger, and direct low pressure refrigerant returning from either the indoor unit heat exchanger or outdoor unit heat exchanger to the inlet of the compressor. In the case where high pressure refrigerant is directed to the outdoor unit heat exchanger, the outdoor unit heat exchanger acts as a condenser and the indoor unit acts as an evaporator, wherein the system rejects heat from the zone to the ambient environment, which is operationally referred to as “cooling mode.” Conversely, in the case where the 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 ambient environment and pumping this heat into the zone, which is operationally referred to as “heating mode.”
In some embodiments, the VCS 1400 can be modeled as a hybrid nonlinear system and the commands that are sent to the actuators are computed using a predictive controller, such as the mixedinteger nonlinear model predictive controller (MINMPC). For example, in some embodiments, the commands that are sent to the valves and/or the fans can only take a discrete set of values, and therefore resulting into a set of binary or integer control input variables for each stage within the mixedinteger control horizon.
In some embodiments, the predictive controller determines an input to the vapor compression system based on the mixedinteger optimal control solution, wherein the input to the vapor compression system includes one or a combination of an indoor unit fan speed, an outdoor unit fan speed, a compressor rotational speed, an expansion valve position, and a flow reversing valve position, and the discrete optimization variables are used to model one or a combination of discrete control decisions, switching in the system dynamics, and integer values for the commands that are sent to the valves and/or to the fans.
In some embodiments, the dynamic behavior of the VCS 1400 can change rapidly or even switch at certain time instances, depending 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 mixedinteger control horizon.
The abovedescribed embodiments of the present invention can be implemented in any of numerous ways. For example, the embodiments may be implemented using hardware, software or a combination thereof. When implemented in software, the software code can be executed on any suitable processor or collection of processors, whether provided in a single computer or distributed among multiple computers. Such processors may be implemented as integrated circuits, with one or more processors in an integrated circuit component. Though, a processor may be implemented using circuitry in any suitable format.
Also, the various methods or processes outlined herein may be coded as software that is executable on one or more processors that employ any one of a variety of operating systems or platforms. Additionally, such software may be written using any of a number of suitable programming languages and/or programming or scripting tools, and also may be compiled as executable machine language code or intermediate code that is executed on a framework or virtual machine. Typically, the functionality of the program modules may be combined or distributed as desired in various embodiments.
Also, the embodiments of the invention may be embodied as a method, of which an example has been provided. The acts performed as part of the method may be ordered in any suitable way. Accordingly, embodiments may be constructed in which acts are performed in an order different than illustrated, which may include performing some acts concurrently, even though shown as sequential acts in illustrative embodiments.
Although the invention has been described by way of examples of preferred embodiments, it is to be understood that various other adaptations and modifications can be made within the spirit and scope of the invention. Therefore, it is the object of the appended claims to cover all such variations and modifications as come within the true spirit and scope of the invention.
Claims
1. A predictive feedback controller for controlling a hybrid dynamical system with nonlinear dynamics and continuous and discrete elements of operation, the predictive feedback controller comprising: at least one processor; and a memory having instructions stored thereon that, when executed by the at least one processor, cause the predictive feedback controller to:
 accept feedback signal including measurements indicative of a current state of the hybrid dynamical system including one or a combination of a current state of the predictive controller, a current state of one or multiple actuators of the hybrid dynamical system, and a current state of outputs of the hybrid dynamical system;
 formulate a mixedinteger nonlinear programming (MINLP) problem optimizing an objective function subject to one or multiple constraints with a solution indicative of a control command for changing the current state of the hybrid dynamical system according to a control objective, wherein the constraints include equality constraints, inequality constraints, or both, and wherein the constraints and the control objective of the MINLP problem include one or multiple nonlinear functions of continuous optimization variables representing the continuous elements of the operation of the hybrid dynamical system and one or multiple linear functions of integer optimization variables representing the discrete elements of the operation of the hybrid dynamical system, such that the MINLP problem is formulated into a separable format ensuring that the discrete elements of the operation are present only in the linear functions of the MINLP problem;
 solve the MINLP problem over multiple iterations of a sequential convexificationbased optimization procedure until a termination condition is met, wherein, to perform an iteration, the predictive feedback controller is configured to perform a partial convexification of a portion of a space of the solution including a current solution guess, wherein the partial convexification produces a convex approximation of the nonlinear functions of the MINLP without approximating the linear functions of the MINLP to produce a partially convexified MINLP; and update the current solution guess by solving a mixedinteger convex programming (MICP) formulation of the partially convexified MINLP problem; and
 submit the control command generated according to the solution of the MINLP problem to the hybrid dynamical system thereby causing a change of the current state of the hybrid dynamical system.
2. The predictive feedback controller of claim 1, wherein the hybrid dynamical system includes a device controlled by manipulation input signals representing the control command, wherein the input signals specify values of the continuous elements of the operation of the hybrid dynamical system and values of the discrete elements of the operation of the hybrid dynamical system, wherein the continuous elements include one or a combination of voltages, pressures, forces, torques, steering angles, velocities, and temperatures, and wherein the discrete elements include one or a combination of energy levels, quantized valve inputs, gear shifts, on/off actuation, lane selection, and obstacle avoidance decision variables.
3. The predictive feedback controller of claim 1, wherein the processor is configured to
 transform the MINLP problem from an original format into the separable format by adding one or a combination of additional integer optimization variables and additional continuous optimization variables to the original format, thereby ensuring that the discrete elements of the operation are present only in the linear functions of the MINLP in separable format.
4. The predictive feedback controller of claim 1, wherein to perform the partial convexification, the processor is configured to:
 identify if a function in the objective function or the constraints of the MINLP problem is a function of the integer optimization variables or the continuous optimization variables;
 compute a convex approximation of the smooth nonlinear function in a local neighborhood of the current solution guess within the solution space of the MINLP problem if the function is a smooth nonlinear function of the continuous optimization variables; and
 avoid the convex approximation by preserving the integer variables, lying within a nonconvex portion of the solution space of the MINLP, if the function is a linear function of integer optimization variables.
5. The predictive feedback controller of claim 4, wherein the convex approximation is computed by a local constraint linearization or a linearquadratic objective approximation of the smooth nonlinear function in a local neighborhood of the current solution guess within the solution space of the MINLP problem, by evaluating one or multiple first and/or higherorder directional derivatives for the smooth nonlinear function, using symbolic differentiation, numerical differentiation or algorithmic differentiation.
6. The predictive feedback controller of claim 1, wherein to update the current solution guess of the iteration, the predictive feedback controller is configured to:
 solve the MICP problem of the partially convexified MINLP to compute a search direction for the continuous optimization variables and a search direction for the integer optimization variables;
 update current values of the integer optimization variables toward the search direction for the integer optimization variables with a step satisfying constraints on each of the integer optimization variables being an integervalued; and
 update current values of the continuous optimization variables toward the search direction for the continuous optimization variables with a step size value governed by feasibility and optimality conditions of the MINLP problem.
7. The predictive feedback controller of claim 6, wherein the step size value is between zero and one and is selected for the iteration based on a merit function balancing optimality and feasibility for the solution of the MINLP problem, such that a value for the merit function decreases between at least two iterations of the sequential convexificationbased optimization procedure.
8. The predictive feedback controller of claim 7, wherein the step size value is selected based on an iterative line search procedure, starting from an initial step size value of one, in which the predictive feedback controller performs one or multiple iterations of the iterative line search procedure to select the step size value and each iteration is configured to:
 test whether a decrease condition is satisfied for a merit function balancing optimality and feasibility for the solution of the MINLP problem;
 terminate the iterative line search procedure and use the step size value to update current values of the continuous optimization variables toward the search direction for the continuous optimization variables, if the sufficient decrease condition is satisfied; and otherwise
 decrease the step size value and continue the iterative line search procedure is continued.
9. The predictive feedback controller of claim 7, wherein the merit function is a combination of the objective function and one or multiple penalty functions, which are applied to a violation for each of the equality constraints and a violation for each of the inequality constraints in the MINLP problem.
10. The predictive feedback controller of claim 6, wherein the MICP problem of the partially convexified MINLP includes one or multiple trustregion inequality constraints, ensuring that the search direction for the continuous optimization variables is below a trustregion radius value.
11. The predictive feedback controller of claim 10, wherein to update the current solution guess of the iteration, the processor is configured to:
 compute a ratio of an actual to a predicted reduction for a value of a merit function balancing optimality and feasibility for the MICP solution of the partially convexified MINLP problem; and
 update current values of the integer and continuous optimization variables toward the search direction for the integer and continuous optimization variables if the ratio of an actual to a predicted reduction for a value of the merit function is above a threshold.
12. The predictive feedback controller of claim 11, wherein to update the trustregion radius value based on the solution for the MICP problem of the partially convexified MINLP, the processor is configured to:
 decrease the trustregion radius value, ensuring that the search direction for the continuous optimization variables is reduced in one or multiple of the next iterations if the ratio of an actual to a predicted reduction for a value of the merit function or the search direction for the continuous optimization variables is below a threshold; and otherwise
 increase the trustregion radius value, ensuring that the search direction for the continuous optimization variables is increased in one or multiple of the next iterations.
13. The predictive feedback controller of claim 1, wherein a branchandbound (B&B) optimization method is used to compute a globally optimal solution for the MICP problem of the partially convexified MINLP.
14. The predictive feedback controller of claim 13, wherein the MICP problem is a mixedinteger linear programming (MILP) subproblem, a mixedinteger quadratic programming (MIQP) subproblem, a mixedinteger quadratically constrained quadratic programming (MIQCQP) subproblem or a mixedinteger secondorder cone programming (MISOCP) subproblem.
15. The predictive feedback controller of claim 6, wherein the current values of the integer optimization variables are fixed after one or multiple iterations of the sequential convexificationbased optimization procedure, and to update the current solution guess in one or multiple iterations of the sequential convexificationbased optimization procedure, the predictive feedback controller is configured to:
 solve a convex programming (CP) problem of the partially convexified MINLP to compute a search direction for the continuous optimization variables, given the fixed current values of the integer optimization variables; and
 update current values of the continuous optimization variables toward the search direction for the continuous optimization variables with a step size value governed by feasibility and optimality conditions of the MINLP problem.
16. The predictive feedback controller of claim 15, wherein the CP problem is a convex linear programming (LP), a convex quadratic programming (QP), convex quadratically constrained quadratic programming (QCQP) or a convex second order cone programming (SOCP) subproblem.
17. The predictive feedback controller of claim 1, wherein a homotopytype penalty method is used to penalize changes in a value for one or multiple of the integer optimization variables, based on a positive weight value, between at least two iterations of the sequential convexificationbased optimization procedure.
18. The predictive feedback controller of claim 17, wherein the homotopytype penalty method adds one or multiple penalty terms in the objective function of the MICP problem of the partially convexified MINLP, in order to accelerate convergence of the sequential convexificationbased optimization procedure to a feasible and/or optimal solution of the MINLP problem.
19. The predictive feedback controller of claim 17, wherein the homotopytype penalty method adds one or a combination of multiple linear and/or smooth nonlinear inequality constraints to the inequality constraints of the MICP problem of the partially convexified MINLP, in order to accelerate convergence of the sequential convexificationbased optimization procedure to a feasible and/or optimal solution of the MINLP problem.
20. The predictive feedback controller of claim 1, wherein the at least one processor is at least one parallel processor and the solution to the MICP problem of the partially convexified MINLP can be found by searching for a globally optimal solution of the MICP problem in different regions of the solution space in parallel.
21. The hybrid dynamical system including the predictive feedback controller of claim 1, wherein the control command generated by the predictive feedback controller specifies a target state of the hybrid dynamical system, the hybrid dynamical system comprises:
 a tracking controller configured to generate one or multiple control inputs to the actuators of the hybrid dynamical system to reduce the error between the current state of the hybrid dynamical system and the target state of the hybrid dynamical system.
22. The hybrid dynamical system of claim 21, wherein the tracking controller is a PID controller or a model predictive controller (MPC).
23. The predictive feedback controller of claim 1, wherein the predictive feedback controller is implemented using a mixedinteger nonlinear model predictive control (MINMPC), wherein the MINMPC computes the control signal based on current state of the system and control command, and wherein the MINMPC computes a control solution that comprises a sequence of future optimal discrete and continuous control inputs over a prediction time horizon of the hybrid dynamical system, by solving a constrained mixedinteger nonlinear optimization problem at each control time step.
24. The predictive feedback controller of claim 23, wherein the solution of the constrained mixedinteger nonlinear optimization problem at a control time step is used as an initial solution guess for the sequential convexificationbased optimization procedure to compute a solution to the constrained mixedinteger nonlinear optimization problem at the next control time step.
25. The predictive feedback controller of claim 1, wherein the system is a vehicle, and wherein the predictive feedback controller determines an input to the vehicle based on the MINLP solution, wherein the input to the vehicle includes one or a combination of an acceleration of the vehicle, an engine torque of the vehicle, brake torques, and a steering angle, and the discrete optimization variables are used to model one or a combination of discrete control decisions, switching in the system dynamics, gear shifting, and obstacle avoidance constraints.
26. The predictive feedback controller of claim 1, wherein the system is a spacecraft, and wherein the predictive feedback controller determines an input to the spacecraft based on the MINLP solution, wherein the input to the spacecraft actuates one or a combination of thrusters and momentum exchange devices, and the discrete optimization variables are used to model one or a combination of discrete control decisions, switching in the system dynamics, integer values for the thruster commands, and obstacle avoidance constraints.
27. The predictive feedback controller of claim 1, wherein the system is a vapor compression system, and wherein the predictive feedback controller determines an input to the vapor compression system based on the MINLP solution, wherein the input to the vapor compression system includes one or a combination of an indoor unit fan speed, an outdoor unit fan speed, a compressor rotational speed, an expansion valve position, and a flow reversing valve position, and the discrete optimization variables are used to model one or a combination of discrete control decisions, switching in the system dynamics, and integer values for the commands that are sent to the valves and/or to the fans.
Type: Application
Filed: Nov 18, 2021
Publication Date: May 18, 2023
Applicant: Mitsubishi Electric Research Laboratories, Inc. (Cambridge, MA)
Inventors: Rien Quirynen (Somerville, MA), Stefano Di Cairano (Newton, MA)
Application Number: 17/455,472