US20190107072A1 - Hybrid partial and full step quadratic solver for model predictive control of diesel engine air path flow and methods of use - Google Patents

Abstract

Methods and systems for use of model predictive control (MPC) controllers utilizing hybrid, quadratic solvers to solve a linear feasibility problem corresponding to a nonlinear problem for an internal combustion engine plant such as a diesel engine air path. The MPC solves a convex, quadratic cost function having optimization variables and constraints and directs the plant per the output solutions to optimize plant operation while adhering to regulations and constraints. The problem includes a combination of iterative and direct calculations in the primal space depending on whether a partial step (iterative) or a full step (direct) is attempted. Further, primal and dual space array matrices are pre-computed and stored offline and are retrieved via use of a unique identifier associated with a specific active set for a set of constraints. Such hybrid and/or offline calculations allow for a reduction in computational power while still maintaining accuracy of solution results.

Description

CROSS REFERENCE TO RELATED APPLICATIONS
• The present application is a divisional of and claims the benefit of U.S. Non-Provisional patent application Ser. No. 15/185,664, filed Jun. 17, 2016, and entitled “Hybrid Partial and Full Step Quadratic Solver For Model Predictive Control of Diesel Engine Air Path Flow and Methods of Use,” the entirety of which is incorporated by reference.
TECHNICAL FIELD
• The present specification generally relates to application of quadratic solvers to solve a linear feasibility problem corresponding to a nonlinear programming formulation and, more specifically, to application of hybrid, e.g., partial and full step, quadratic solvers to solve a linear feasibility problem corresponding to a nonlinear problem of an internal combustion engine plant (e.g., a diesel engine air path).
BACKGROUND
• In internal combustion engines, an amount of air supplied to engine cylinders may be manipulated by engine components. For example, in modern diesel engines, variable geometry turbines (VGT) may be used to increase an amount of air supplied to engine cylinders by varying an angle of turbine stator inlet vanes such that the amount of supplied air is changed.
• Such modern diesel engines typically balance providing optimum performance and fuel economy while meeting stringent federal regulations on emissions, such as constraints on particulate matter and nitrogen oxides. To meet these requirements, many diesel engines having a VGT also use an exhaust gas recirculation (EGR) valve having a variable controlled position. The EGR valve re-circulates varying amounts of engine exhaust gases back into the engine cylinders to allow for both a more complete combustion and reduced engine emissions.
• Such engines operate over a large range of operating conditions, which may include, for example, engine speed, fuel usage, and engine load, among other conditions, and one or more controllers are embedded in an engine control unit (ECU) to control various engine actuators in response to sensors that detect engine performance. The ECU works to optimize engine performance, emissions, and other such outputs via use of, for example, quadratic solvers.
• Accordingly, a need exists for alternative quadratic solvers to provide reduced computational time and an increased accuracy and methods of use of such quadratic solvers.
SUMMARY
• In one embodiment, a method for controlling an internal combustion engine having a variable geometry turbine (VGT), an exhaust gas recirculation (EGR) valve, and an EGR throttle may include solving a linear quadratic problem with a predictive model comprising an updating algorithm in order to determine: (i) a requested optimized VGT lift that meets one or more constraints; and (ii) a requested optimized EGR valve flow rate that meets the one or more constraints, generating the requested optimized VGT lift responsive to an engine intake manifold pressure by controlling the VGT, and generating the requested optimized EGR valve flow rate responsive to an EGR rate by controlling the EGR valve and the EGR throttle. The solving the linear quadratic problem may include determining whether to take a determined step comprising one of a primal partial step and a primal full step at each iteration, and taking the determined step at each iteration until the linear quadratic problem is solved by the updating algorithm. Taking the primal partial step may include performing an iterative calculation, and taking the primal full step may include performing a direct calculation.
• In another embodiment, a method for controlling an internal combustion engine having, within an air path of the engine, a variable geometry turbine (VGT), an exhaust gas recirculation (EGR) valve, and an EGR throttle may include formulating a constrained optimization problem for a model predictive control (MPC) controller controlling the air path based on a linear model, one or more constraints, and associated dual space and primal space matrix arrays, the linear model comprising a convex, quadratic, time-varying cost function in dual and primal space, and each array associated with a unique active set list comprising a first combination of the one or more constraints. The method may further include solving the constrained optimization problem to determine a solution, updating the constrained optimization problem with an updated active set list, repeating the solving and formulating steps until all possible active set lists of the one or more constraints are satisfied to generate a requested optimized VGT lift and a requested optimized EGR valve flow rate, each of which meets the one or more constraints, to control the air path, and implementing the solution with respect to the air path. The requested optimized VGT lift may be generated responsive to an engine intake manifold pressure by controlling the VGT, and the requested optimized EGR valve flow rate may be generated responsive to an EGR rate by controlling the EGR valve and the EGR throttle.
• In yet another embodiment, a system for controlling an internal combustion engine having, within an air path of the engine, a variable geometry turbine (VGT), an exhaust gas recirculation (EGR) valve, and an EGR throttle may include a processor communicatively coupled to a non-transitory computer storage medium, wherein the non-transitory computer storage medium stores instructions that, when executed by the processor, cause the processor to formulate a constrained optimization problem for a model predictive control (MPC) controller controlling the air path based on a linear model, one or more constraints, and associated dual space and primal space matrix arrays, the linear model comprising a convex, quadratic, time-varying cost function in dual and primal space, and each array associated with a unique active set list comprising a first combination of the one or more constraints. The non-transitory computer storage medium may store further instructions that, when executed by the processor, cause the processor to solve the constrained optimization problem to determine a solution, update the constrained optimization problem with an updated active set list, repeat the solving and formulating steps until all possible active set lists of the one or more constraints are satisfied to generate a requested optimized VGT lift and a requested optimized EGR valve flow rate, each of which meets the one or more constraints, to control the air path, and implement the solution with respect to the air path,. The requested optimized VGT lift may be generated responsive to an engine intake manifold pressure by controlling the VGT, and the requested optimized EGR valve flow rate may be generated responsive to an EGR rate by controlling the EGR valve and the EGR throttle.
• These and additional features provided by the embodiments described herein will be more fully understood in view of the following detailed description, in conjunction with the drawings.
BRIEF DESCRIPTION OF THE DRAWINGS
• The embodiments set forth in the drawings are illustrative and exemplary in nature and not intended to limit the subject matter defined by the claims. The following detailed description of the illustrative embodiments can be understood when read in conjunction with the following drawings, where like structure is indicated with like reference numerals and in which:
• FIG. 1 is a schematic illustration of a diesel engine controlled by a model predictive controller, according to one or more embodiments shown and described herein;
• FIG. 2 is a schematic block diagram showing inputs and outputs to the model predictive controller, according to one or more embodiments shown and described herein;
• FIG. 3 is a schematic block diagram of the engine ofFIG. 1 using the model predictive controller ofFIG. 2, according to one or more embodiments shown and described herein;
• FIG. 4 is an illustration of a flow diagram depicting offline and online sequence steps and operations of a model predictive control method, according to one or more embodiments shown and described herein;
• FIG. 5 is an illustration of a flow diagram of a hybrid quadratic solver algorithm applied with respect to the model predictive control method ofFIG. 4, according to one or more embodiments shown and described herein;
• FIG. 6A is a graphical illustration of a convex cost function problem in a three-dimensional primal space, according to one or more embodiments shown and described herein;
• FIG. 6B is a graphical illustration of the convex cost function problem ofFIG. 6A in a two-dimensional primal space and having a primal feasible region, according to one or more embodiments shown and described herein;
• FIG. 7A is a graphical illustration of a first set of iterations in a primal space in an example of a solution to the convex problem ofFIG. 6A that utilizes the flow diagram ofFIG. 5, according to one or more embodiments shown and described herein;
• FIG. 7B is a graphical illustration of a second set of iterations in the primal space following the first set of iterations ofFIG. 7A, according to one or more embodiments shown and described herein;
• FIG. 7C is a graphical illustration of a third set of iterations in the primal space following the second set of iterations ofFIG. 7B, according to one or more embodiments shown and described herein;
• FIG. 8A is a graphical illustration of another primal space example of a solution to the convex problem ofFIG. 6A, which includes a primal problem having four constraints and two optimization variables, according to one or more embodiments shown and described herein;
• FIG. 8B is a graphical illustration of a first step in solving the primal problem ofFIG. 8A, according to one or more embodiments shown and described herein;
• FIG. 8C is a graphical illustration of a second step in solving the primal problem ofFIG. 8B, according to one or more embodiments shown and described herein;
• FIG. 8D is a graphical illustration of a third step in solving the primal problem ofFIG. 8C, according to one or more embodiments shown and described herein;
• FIG. 9 illustrates a performance of engine speed controlled by a MPC controller applying the method ofFIGS. 4-5, according to one or more embodiments shown and described herein; and
• FIG. 10 illustrates a maximum error output comparison of the data ofFIG. 9 against transient bench data, according to one or more embodiments shown and described herein.
DETAILED DESCRIPTION
• Referring generally to the figures, embodiments of the present disclosure are directed to use of MPC controllers utilizing hybrid, e.g., partial and full step, quadratic solvers to solve a linear feasibility problem corresponding to a nonlinear problem for an internal combustion engine plant (e.g., the plant being a diesel engine air path) that has engine operating parameters as constraints. Generally, the MPC methods described herein solve a convex, quadratic cost function having optimization variables and constraints and direct operation of the plant per output solutions to optimize plant operation while adhering to regulations and constraints. The problem is time-varying, though is solvable as a discrete time invariant model, and can include a combination of iterative and direct calculations in the primal space depending on whether a partial step (utilizing the iterative calculation) or a full step (utilizing the direct calculation) is attempted. Further, primal and dual space array matrices are pre-computed, stored offline, and are retrieved via use of a unique identifier associated with a specific active set for a set of constraints, as will be described in greater detail further below. Such hybrid and/or offline calculations allow for a reduction in computational power while still maintaining accuracy of solution results for implementation by the plant.
• Various embodiments of the diesel engine air path plant and MPC operations and methods for use of such plant operations via MPC controllers are described in detail herein. It should be understood that the algorithms described herein may be applied to plants other than diesel engine air paths as will be apparent to those of ordinary skill in the art.
• FIG. 1 depicts aninternal combustion engine20, which is, as a non-limiting example, a diesel engine. Theengine20 includes anengine block22 that houses a plurality ofcylinders24. Afuel rail26 is connected to a fuel supply (not shown) and supplies diesel fuel to a plurality offuel injectors28, and eachcylinder24 is provided with afuel injector28.
• Anintake manifold30 is coupled to the plurality ofcylinders24 to supply intake air to eachcylinder24. Also coupled to theintake manifold30 is an intakemanifold pressure sensor32 to measure intake manifold air pressure. Combustion gases are carried away from the plurality ofcylinders24 and theengine block22 by anexhaust manifold34.
• Abypass path29 between theintake manifold30 and theexhaust manifold34 has a coupledEGR valve40 to re-circulate a portion of the exhaust gases from theexhaust manifold34 back into theintake manifold30 for supply to the plurality ofcylinders24. Along with theEGR valve40, anEGR cooler42 may be coupled in thebypass path29. As described above, the EGR valve re-circulates varying amounts of engine exhaust gases back into the engine cylinders to allow for both a more complete combustion and reduced engine emissions. The amount theEGR valve40 is opened controls an amount of engine exhaust gases which are able to re-circulate through thebypass path29 from theexhaust manifold34 back into theintake manifold30. TheEGR cooler42 assists to help prevent theEGR valve40 from overheating, which may otherwise lead to an increased wear and tear.
• AnEGR throttle44 to further assist with controlling gas circulation is mounted in an airflow path from acompressor46 of aVGT48. Anintercooler50 to assist with preventing overheating of theEGR throttle44 may be mounted ahead of theEGR throttle44 in the airflow path for intake air. Thecompressor46 increases a pressure of the incoming air. Further, theVGT48 includes turbine input vanes that may be opened, partially opened, or closed through an angling of the turbine input vanes to control a VGT lift and to allow for the passage of air in through theEGR throttle44 to join with the exhaust gases being re-circulated into theintake manifold30 through thebypass path29. Thus, by controlling an angle of turbine input vanes, theVGT48 controls an intake manifold pressure provided by thecompressor46 of theVGT48. An amount theEGR throttle44 is opened also restricts the amount of air provided through theVGT48 that is able to join with air recirculated through thebypass path29. Also coupled to walls defining thebypass path29 is anEGR rate sensor43 to measure EGR rate (such as a fraction of re-circulated air versus fresh air) as it is affected by theEGR valve40 and/or theEGR throttle44. Another measurement, EGR flow rate or EGR flow, may refer to an amount of mass re-circulated air flow through theEGR valve40.
• In embodiments described herein, a model predictive controller utilizing a MPC for theengine20 uses a plurality of control inputs. The control inputs may be, for example, and as shown inFIG. 2, anintake manifold pressure62 measured or calculated from signals generated by the intakemanifold pressure sensor32 and/or anEGR rate64 measured or calculated from signals generated by theEGR rate sensor43. One or more set point maps may prescribe set-points for the control inputs such as intake pressure and EGR rate. Further, Kalman filter estimates based on measurements of control inputs such as intake pressure, ECU estimated EGR rate, and mass airflow (MAF) may be used.
• In the MPC controller methods described herein, and as described in greater detail further below, a single solution is solved for a convex problem based on dual and primal space arrays, which may be pre-computed and stored offline. As a system that is utilized for the MPC controller methods described herein,FIG. 3 shows anECU70, which controls the plant (e.g., the diesel engine air path) ofFIGS. 1-2. TheECU70 includes aprocessor61 that may execute a computer program that is tangibly embodied on a computer usable medium and includes instructions that are executable by the processor to direct the MPC controller as described below. TheECU70 may including a central processing unit (CPU) that is any type of device or devices capable of manipulating or processing information. The CPU may be practiced with a single or multiple processors.
• TheECU70 may be communicatively coupled to the actuators (e.g., to actuate theEGR valve40, theEGR throttle44, and/or the VGT48),engine20, and sensors (e.g., the intakemanifold pressure sensor32 and/or the EGR rate sensor43) as shown inFIG. 3 through acommunication path65, for example. Thecommunication path65 may be formed from any medium that is capable of transmitting a signal such as, for example, conductive wires, conductive traces, optical waveguides, or the like, or from a combination of mediums capable of transmitting signals. Thecommunication path65 communicatively couples the various components of theengine20 and theECU70. As used herein, the term “communicatively coupled” means that coupled components are capable of exchanging data signals with one another such as, for example, electrical signals via conductive medium, electromagnetic signals via air, optical signals via optical waveguides, and the like.
• As noted above, theECU70 may include one ormore processors61 that can be any device capable of executing machine readable instructions. Accordingly, aprocessor61 may be a controller, an integrated circuit, a microchip, a computer, or any other computing device. Theprocessor61 is communicatively coupled to the other components ofFIG. 3 by thecommunication path65. Accordingly, thecommunication path65 may communicatively couple any number ofprocessors61 with one another, and allow the modules coupled to the communication path to operate in a distributed computing environment. Specifically, each of the modules can operate as a node that may send and/or receive data.
• TheECU70 also includes amemory component63, which is coupled to thecommunication path65 and communicatively coupled to theprocessor61. Thememory component63 may be a non-transitory computer readable medium and may be configured as a nonvolatile computer readable medium. Thememory component63 may comprise RAM, ROM, flash memories, hard drives, or any device capable of storing machine readable instructions such that the machine readable instructions can be accessed and executed by theprocessor61. The machine readable instructions may comprise logic or algorithm(s) written in any programming language such as, for example, machine language that may be directly executed by the processor, or assembly language, object-oriented programming (OOP), scripting languages, microcode, etc., that may be compiled or assembled into machine readable instructions and stored on the memory component. Alternatively, the machine readable instructions may be written in a hardware description language (HDL), such as logic implemented via either a field-programmable gate array (FPGA) configuration or an application-specific integrated circuit (ASIC), or their equivalents. Accordingly, the methods described herein may be implemented in any conventional computer programming language, as pre-programmed hardware elements, or as a combination of hardware and software components. TheECU70 further includes additional storage ordatabases67 to store components such as off-line pre-computed matrices, as described in greater detail further below. Thememory component63 may include machine readable instructions that, when executed by theprocessor61, cause theprocess61 to perform the functions of theECU70, operating as an MPC controller.
• TheECU70 includes thenetwork interface hardware69 for communicatively coupling theECU70 with a computer network71. Thenetwork interface hardware69 is coupled to thecommunication path65 such that thecommunication path65 communicatively couples thenetwork interface hardware69 to other modules of theECU70. Thenetwork interface hardware69 can be any device capable of transmitting and/or receiving data via a wireless network. Accordingly, thenetwork interface hardware69 can include a communication transceiver for sending and/or receiving data according to any wireless communication standard. For example, the network interface hardware can include a chipset (e.g., antenna, processors, machine readable instructions, etc.) to communicate over wired and/or wireless computer networks such as, for example, wireless fidelity (Wi-Fi), WiMax, Bluetooth, IrDA, Wireless USB, Z-Wave, ZigBee, or the like.
• The quadratic programming (QP) solver as described herein utilizes, in a first embodiment, a QPKWIK algorithm with a derivation as set forth herein for controlling a plant of a diesel engine air path. Utilizing the QPKWIK algorithm for controlling a diesel engine air path plant allows for finding an MPC solution for large scale process optimization problems that include many variables and time-varying constraints but only a few degrees of freedom, which helps reduce computational time. The optimization variables as time varying control inputs may be, for example, aVGT lift66 and a EGRvalve flow rate68. Time varying linear terms or variables may be, for example, control inputs such as anintake manifold pressure62 and/or anEGR rate64. The QPKWIK algorithm uniquely only requires that an inverse Cholesky factor of a Hessian matrix to be supplied, which inverse Cholesky factor is obtained at each iteration directly using a factorized inverse BFGS formula (e.g., a quasi-Newton update formula). A resulting quadratic program sub-problem has a resulting solution having degrees of freedom of the second order rather than a third order, requiring less computational time than a higher order problem, for example.
• Below is a derivation of the QPKWIK algorithm.
• The following convex quadratic problem serves as a convex function to solve by theECU70 while satisfying the associated problem constraints:
• $\begin{array}{cc}\underset{}{\min }J\left(\right)=\frac{1}{2}\left(T^{}+{\mathscr{H}}^T\right)s.t.\le L_{}\le \le U_{}& \left(\mathrm{Equation}3.1\right)\end{array}$
• At each time step k,

  

  is a linear term in the cost function (time varying);

  

  is a quadratic term in the cost function (constant);

  

  is a constraint matrix (constant);

  

  is a constraint vector (time varying); and

  

  is one or more optimization variables (e.g., time varying control inputs). Assumptions are made that J is convex and that the constraints are linear. Equation 3.1 can be simplified by taking into account only active (feasible) constraints to the following:
• $\begin{array}{cc}\underset{}{\min }J\left(\right)=\frac{1}{2}\left(T^{}+{\mathscr{H}}^T\right)s.t.\le & \left(\mathrm{Equation}3.2\right)\end{array}$
• Active constraints indicate those constraints that are possible and plays a part in finding an optimal solution such that the inequality active constraint holds as an equality constraint, as described in greater detail further below. Using, for optimality, Karush-Kuhn-Tucker (KKT) conditions, which must hold for a solution to be a minimum, and a vector of Lagrange multipliers, the convex problem and constraints are combined to form the following Lagrangian:
• $\begin{array}{cc}ℒ=\frac{1}{2}T^{}+\mathscr{H}-{\mu }^T\left(-\right)& \left(\mathrm{Equation}3.3\right)\end{array}$
• Taking partial derivatives of Equation 3.3 and setting them equal to zero, optimal values for

  

  * and μ* are obtained as follows and as set forth in Equation 3.4 as a dual problem having two sets of equations to satisfy:
• $\begin{array}{cc}\frac{dℒ}{d}={*}^{}-{\mathrm{\mu }}^{*}+\mathscr{H}=0\frac{dℒ}{d\mu }={}^T{\mu }^{*}+=0& \left(\mathrm{Equation}3.4\right)\end{array}$
• Further, for the dual problem to have dual feasibility, all of the Lagrange multipliers μ must be greater than or equal to zero. Those that are equal to zero are inactive constraints, those that are greater than zero are active constraints playing a role in the solution, and those that are less than zero indicate the problem is infeasible and should be investigated. If the optimal value for

  

  satisfies Equation 3.4 above and the Lagrange multipliers μ are zero for all inactive constraints, then the dual cost function

  

  (

  

  , μ) is equal to the primal cost function J(

  

  ) and a complementary slackness condition is upheld such that the solution to the dual problem is the same as the solution to the primal problem.
• Thus, if there are no active constraints, then Equation 3.4 returns the unconstrained minimum (e.g., the global minimum) of the convex problem. Equation 3.4 can be written into matrix form, as shown in Equation 3.5 below:
• $\begin{array}{cc}\left[\begin{array}{cc}& -\\ -{}^T& 0\end{array}\right]\left[\begin{array}{c}{*}^{}\\ {\mu }^{*}\end{array}\right]=-\left[\begin{array}{c}\mathscr{H}\\ \end{array}\right]& \left(\mathrm{Equation}3.5\right)\end{array}$
• Equation 3.5 then is able to become Equation 3.6, which is the solution for the optimization variable(s)

  

  * and Lagrangian muliplier(s) μ* as set forth below:
• $\begin{array}{cc}\left[\begin{array}{c}{*}^{}\\ {\mu }^{*}\end{array}\right]=\left[\begin{array}{cc}-H& D\\ D^T& -U\end{array}\right]\left[\begin{array}{c}\mathscr{H}\\ \end{array}\right]& \left(\mathrm{Equation}3.6\right)\end{array}$
• where, as set forth in Equation 3.7 below,
• 
  H=

  

  ⁻¹−

  

  ⁻¹

  

  (

  

  ^T

  

  ⁻¹

  

  )⁻¹

  

  ^T

  

  ⁻¹
• 
  D=

  

  ⁻¹

  

  (

  

  ^T

  

  ⁻¹

  

  )⁻¹
• 
  U=−(

  

  ^T

  

  ⁻¹

  

  )⁻¹  (Equation 3.7).
• As

  

  ⁻¹

  

  is repeated many times, Equation 3.7 may be simplified by using (1) a Cholesky decomposition to make

  

  ⁻¹easily invertible and by using (2) QR factorization to allow for inverting a non-square matrix. With respect to the Cholesky decomposition, for a positive definite matrix such as

  

  , the matrix can be decomposed into two triangular matrices L using the following rule:
• 
  

  

  =LL^T  (Equation 3.8)
• Using Cholesky decomposition on

  

  ⁻¹

  

  results in the following Equation 3.9:
• 
  

  

  ⁻¹

  

  =LL^T−1

  

  =L^−TL⁻¹

  

    (Equation 3.9)
• Using QR factorization, which is a decomposition that allows a n×m matrix N to be decomposed into an orthogonal matrix and an upper triangular matrix R, the following rule can be used:
• $\begin{array}{cc}\begin{array}{c}N=\hat{Q}\left[\begin{array}{c}R\\ 0\end{array}\right]\\ =\left[\begin{array}{cc}{\hat{Q}}_1& {\hat{Q}}_2\end{array}\right]\left[\begin{array}{c}R\\ 0\end{array}\right]\\ ={\hat{Q}}_1R+{\hat{Q}}_20\\ ={\hat{Q}}_1R\end{array}& \left(\mathrm{Equation}3.10\right)\end{array}$
• Applying the QR factorization of Equation 3.10 to Equation 3.9 results in the following:
• $\begin{array}{cc}\begin{array}{c}{Q}^{-1}=L^{-T}\left[\begin{array}{cc}{\hat{Q}}_1& {\hat{Q}}_2\end{array}\right]\left[\begin{array}{c}R\\ 0\end{array}\right]\\ =\left[\begin{array}{cc}{L}^{-T}{\hat{Q}}_1& L^{-T}{\hat{Q}}_2\end{array}\right]\left[\begin{array}{c}R\\ 0\end{array}\right]\\ =\left[\begin{array}{cc}{T}^C& T^U\end{array}\right]\left[\begin{array}{c}R\\ 0\end{array}\right]\end{array}& \left(\mathrm{Equation}3.11\right)\end{array}$
• where
• 
  T^C=L^−T{tilde over (Q)}₁
• 
  T^U=L^−T{tilde over (Q)}₂  (Equation 3.12)
• Combining Equations 3.7 and 3.11 lead to matrices, as set forth below, that are used to the solve the KKT condition and are redefined as follows:
• 
  H=T^U(T^U)^T
• 
  D=T^CR^−T
• 
  U=−(R⁻¹R^T)⁻¹  (Equation 3.13)
• Thus, applying Equation 3.13 to Equation 3.6, the solution to the problem is given as follows:
• 
  

  

  =−T^U(T^U)^T

  

  +T^CR^−T

  
• 
  μ*=R⁻¹(T^C)^T

  

  +R⁻¹R^−T

  

    (Equation 3.14A)
• Because the algorithm solves the KKT problem based on a current active set (denoted as subscript L below), Equation 3.14A can be rewritten as Equation 3.14B below:
• 
  

  

  _L=−T_L^U(T_L^U)^T

  

  +T_L^CR_L^−T

  
• 
  μ_L=R_L⁻¹(T_L^C)^T

  

  +R_L⁻¹R_L^−T

  

    (Equation 3.14B)
• As solving the above problem may be computationally costly to do for each iteration, the algorithm may be simplified further and may calculate how to move from an active set L to a new active set L+1, as follows:
• 
  

  

  _L+1=

  

  _L+T_L^U(T_L^U)^TV_knextt
• 
  μ_L+1+μ_L−R_L⁻¹(T_L^C)^Tv_knextt  (Equation 3.15)
• In the above Equation 3.15, t is a minimum length the algorithm can take and still maintain dual feasibility (and is calculated in both primal and dual spaces), and v_knextis the next constraint to be added to a current constraint set list. Additionally, where z is an array that is indicative of a search direction in primal space, and r is an array that is indicative of a search direction in dual space, the following may be used:
• 
  z=T_L^U(T_L^U)T
• 
  r=R_L⁻¹(T_L^C)^T  (Equation 3.16)
• The final QPKWIK equations may then be, in primal and dual space respectively, simplified to the following:
• 
  

  

  _L+1=

  

  _L+zv_knextt
• 
  μ_L+1=μ_L−rv_knextt  (Equation 3.17)
• A minimum value of t (which is calculated for both dual space and primal space) is chosen to ensure the solution stays within a feasible region.
• In a second embodiment, and as described herein, an HQPKWIK algorithm is utilized as a quadratic solver for model predictive control (MPC) of diesel engine air path flow. The HQPKWIK algorithm uses a hybrid model that differs from the QPKWIK algorithm, as will be described in greater detail below. First, the HQPKWIK algorithm utilizes pre-computed primal and dual space arrays that are computed offline (and a lookup algorithm to retrieve the necessary arrays), whereas the QPKWIK algorithm calculates such arrays online such that the iterative, online calculations (as shown by the derivation below) utilize more computational power and read-only memory (ROM) than the HQPKWIK algorithm.
• For example, to determine which matrix to extract from a set of pre-computed arrays, the current set list L is converted to an identifier using the following pseudo-code as a lookup algorithm:

• 	for i = 1 ... l

  if Ll>0

  identifier = sum(identifier) + 2l−i

  end

  	end
  	return identifier

• Once a unique identifier is determined for a given set list L, a pair of z and r (primal and dual space) matrices relating to the set list and associated with the unique identifier may be retrieved. As a non-limiting example, for a given problem with 5 constraints such that 1=5, a current active set list L may haveactive constraints 1, 3, and 5 such that L={1,0,3,0,5}. Using the lookup algorithm, the identifier is able to be determined as follows:
• 
  identifier=2⁵⁻¹+2⁵⁻³+2⁵⁻⁵=2⁴+2²+2°=16+4+1=21
• Thus, the unique identifier is determined to be adecimal number21 from a summation of converted binary numbers to result in a decimal value, and the primal z and dual r space matrices corresponding to this unique identifier may be retrieved. Alternatively, as set forth further below in Example 2, if a constraint is active, it is able to be assigned a value of 1 from which a binary number may be determined and converted to a decimal number to find an associated unique identifier.
• Second, the HQPKWIK algorithm modifies the updating QPKWIK algorithm to improve numerical solutions by incorporating a partial and full step hybrid model utilizing both iterative and direct next step calculation approaches for the primal space. The QPKWIK algorithm, by contrast, utilizes iterative next step calculations for the primal space that may result at times in out of bound errors and infeasible solutions.
• For example, in the QPKWIK algorithm, only iterative computations in primal space are used (as shown by use of v_knext) through the following primal space solution from Equation 3.15:
• 
  

  

  _L+1=

  

  _L+T_L^U(T_L^U)^Tv_knexttPrimal Full Step
• This may result in a computational expense plus possible errors for full step solutions that are greater than an acceptable error of 1e-6 and are outside of the boundaries for a feasible solution. Thus, the hybrid QPKWIK solution, or HQPKWIK, allows for partial steps (also referable to as half steps, though not necessarily having to be exact half portions of a full step) to be taken rather than a full primal step. For example, the partial steps involve iterative calculations (see Equation 3.15) and the primal full step involves direct calculations (as shown by Equation 3.14B) in which the equation is solved directly and doesn't assume where the step is going to be. Thus, Equation 3.15 for primal space is modified for the HQPKWIK algorithm to be as follows:
• $\begin{array}{cc}{u}_{L+1}=\{\begin{array}{cc}{u}_L+{T_L^U\left(T_L^U\right)}^Tv_{\mathrm{knext}}t& \mathrm{Primal}\mathrm{Partial}\left(\mathrm{Half}\right)\mathrm{Step}\\ \begin{array}{c}-{T_{L+1}^U\left(T_{L+1}^U\right)}^T\mathscr{H}+\\ T_{L+1}^CR_{L+1}^{-T}{}_{L+1}\end{array}& \mathrm{Primal}\mathrm{Full}\mathrm{Step}\end{array}\begin{array}{cc}{\mu }_{L+1}={\mu }_L-{R_L^{-1}\left(T_L^C\right)}^Tv_{\mathrm{knext}}t& \mathrm{Dual}\mathrm{Step}\end{array}& \left(\mathrm{Equation}3.18A\right)\end{array}$
• Using the definitions of primal and dual space matrices z and r, respectively, this simplifies to the following equation:
• $\begin{array}{cc}{}_{L+1}=\{\begin{array}{cc}{}_L+{\mathrm{zv}}_{\mathrm{knext}}t& \mathrm{Primal}\mathrm{Partial}\left(\mathrm{Half}\right)\mathrm{Step}\\ \begin{array}{c}-{T_{L+1}^U\left(T_{L+1}^U\right)}^T\mathscr{H}+\\ T_{L+1}^CR_{L+1}^{-T}{}_{L+1}\end{array}& \mathrm{Primal}\mathrm{Full}\mathrm{Step}\end{array}\begin{array}{cc}{\mu }_{L+1}={\mu }_L-{\mathrm{rv}}_{\mathrm{knext}}t& \mathrm{Dual}\mathrm{Step}\end{array}& \left(\mathrm{Equation}3.18B\right)\end{array}$
• Which is further able to be simplified to the following:
• $\begin{array}{cc}{}_{L+1}=\{\begin{array}{cc}{}_L+{\mathrm{zv}}_{\mathrm{knext}}t& \mathrm{Primal}\mathrm{Partial}\left(\mathrm{Half}\right)\mathrm{Step}\\ -z\mathscr{H}+r^T{}_{L+1}& \mathrm{Primal}\mathrm{Full}\mathrm{Step}\end{array}\begin{array}{cc}{\mu }_{L+1}={\mu }_L-{\mathrm{rv}}_{\mathrm{knext}}t& \mathrm{Dual}\mathrm{Step}\end{array}& \left(\mathrm{Equation}3.18C\right)\end{array}$
• FIG. 4 illustrates a flow chart of anoperation400 of an MPC model utilizing the offline and online aspects of the HQPKWIK algorithm, andFIG. 5 illustrates the specific HQPKWIK algorithm that is involved in a solving a constrained optimization problem ofstep416 of the flow chart ofFIG. 4.
• Referring toFIG. 4, theoperation400 of the MPC model utilizing the HQPKWIK algorithm has anoffline portion402 including steps404-410. Instep404, theoperation400 develops and calibrates a nonlinear model of a diesel air path of theengine20. For example, inputs may include VGT position (e.g., theVGT position80 ofFIGS. 2-3) and EGT valve angle (e.g., theEGR valve position82A ofFIGS. 2-3), and outputs may respectively include exhaust/intake pressure62 andEGT rate64 ofFIGS. 2-3, as well as MAF. The inputs and outputs may be generalized and not necessarily restricted to being rate-based, and a derivative operations to inputs and outputs are thus not required for generalized, non-rate based inputs and outputs.
• Instep406, the nonlinear plant is linearized around an operating point such that a linear model is identified from the nonlinear plant. Instep408, constraints on the system are identified, examples of which may include but are not limited to a maximum MAF overshoot, a maximum EGR rate, a maximum/minimum EGR valve command (such as a valve open command) indicating an engine operating range, a maximum/minimum VGT lift closed command, and/or a maximum/minimum EGR flow command. Instep410, offline matrices in the dual and primal search direction are pre-computed and stored. These matrices are used to solve a sub-problem created by using the various combinations of constraints that follow the HQPKWIK algorithm, as described in greater detail below, and by applying the equations of the HQPKWIK algorithm that stem from the QPKWIK derivation as described above.
• Theoperation400 of the MPC model utilizing the HQPKWIK algorithm further has anonline portion412 that includes steps414-418. Instep414, the MPC problem (e.g., Equation 3.2) is formulated. Example inputs are current constraints, observer estimates, sensor inputs, targets, and the like. Example outputs are matrices that define the optimization problem.
• Instep416, the MPC/constrained optimization problem is solved and a desired control is output. For example, and by application of the HQPKWIK algorithm, Equation 3.2 is solved for a current optimization problem. The HQPKWIK algorithm applied in thisstep416 is described in greater detail with respect toFIG. 5 below. Referring again toFIG. 4, instep418, the desired control outputs are provided to the plant of the diesel air path such that the VGT and EGR valves are controlled via the solution provided instep416. Theoperation400 loops back to thestep414 to repeat the steps414-418 to solve and implement solutions for new, ongoing, and next optimization problem sets.
• Referring toFIG. 5, the HQPKWIK algorithm asreference500 is applied to solve the MPC/constrained optimization problem ofstep416. TheHQPKWIK algorithm500 starts atstep501 and continues on to set the optimization variables to an unconstrained minimum (e.g., a global minimum

  

  ₀*) atstep502.
• For example, the global minimum for the convex quadratic function of
• $J\left(\right)=\frac{1}{2}{}^T\mathrm{}+{\mathscr{H}}^T$
• is found by using the gradient, as follows:
• 
  ∇J(

  

  *)=0⇒

  

  

  *+

  

  =0⇒

  

  *=−

  

  ⁻¹

  
• For example, below is an example problem and solution to minimize a quadratic cost function J(

  

  ) having two optimization variables

  

  ₀₁and

  

  ₀₂:
• Problem:
• $\underset{}{\min }J\left(\right)=\frac{1}{2}\left[\begin{array}{cc}{}_{01}& {}_{02}\end{array}\right]\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]\left[\begin{array}{c}{}_{01}\\ {}_{02}\end{array}\right]+\left[\begin{array}{cc}1& 1\end{array}\right]\left[\begin{array}{c}{}_{01}\\ {}_{02}\end{array}\right]$
• Solution:
• $\nabla J\left({}^{*}\right)=-Q^{-1}\mathscr{H}=-{\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]}^{-1}\left[\begin{array}{c}1\\ 1\end{array}\right]\Rightarrow {}^{*}=\left[\begin{array}{c}{}_{01}^{*}\\ {}_{02}^{*}\end{array}\right]=\left[\begin{array}{c}-1\\ -1\end{array}\right]$
• FIG. 6A shows an example of a graphical representation of the convex quadratic function J(

  

  ) in a three-dimensional primal space, andFIG. 6B shows a graphical representation of the convex quadratic function in a primal space including the global minimum

  

  ₀* (which is the minimum, bottom-most point of the three-dimensional convex cost function ofFIG. 6A).FIG. 6B further shows, as will be described in greater detail further below, an active (possible) set list L={1,2} having two (active) constraints.FIG. 6B shows a first iteration

  

  ₁* alongconstraint1, a second iteration

  

  ₂* alongconstraint2, and a point

  

  _1,2when a set list L={1,2} is an active, possible one. The region above point

  

  _1,2and contained by the set list L={1,2} is the primal feasible region, and point

  

  _1,2* is a minimum feasible solution which may be found via theHQPKWIK algorithm500 as described in greater detail below. Primal feasibility as described herein refers to whether a point lies within a region that satisfies all constraints (and is termed the “primal feasible region”). If the point does lie within the primal feasible region, then the point is feasible, and if it does not, then it is infeasible.
• Further, with respect to active constraints versus inactive constraints as described herein, if an inequality constraint is considered active, then the inequality constraint (e.g., h_i(

  

  )=

  

  _i

  

  +

  

  _i≤0) plays a part in finding the optimal solution

  

  * and holds as an equality constraint, e.g., h_i(

  

  )=0. However, while a constraint may be inactive (as the optimal point

  

  * does not lie on an intersection stemming from the constraint), it may still create a feasible region on which the optimal point

  

  * lies.
• Referring again toFIG. 5, after the unconstrained (global) minimum

  

  ₀* is found instep502, as described above, the most violated (possible) constraint is determined instep504. Instep506, if all the possible constraints are satisfied, the algorithm exits atstep508 to proceed onto the plant (e.g., diesel air path) implementation instep418 ofFIG. 4 with the solution that includes the desired control outputs. If, however, all the possible constraints are not satisfied instep506, the algorithm continues on to step508 to convert an active set list to a unique identifier with respect to a lookup algorithm, as described above, for example. Once the unique identifier is obtained, then the identifier is able to be used to locate specific pre-computed dual and primal space arrays (r and z, respectively, of Equations 3.16-3.17) to be used instep512 to calculate primal and dual step directions and determine step length. Instep514, the Lagrange multipliers and the optimization variables are updated. Instep516, the algorithm determines what type of step is to be taken.
• For example, a primal full step is taken for feasible step directions that are within an error boundary and satisfy dual feasibility. As a result of this step, an appropriate constraint is added instep520 and the algorithm returns to step504 to determine the next most violated constraint and to repeat the steps as described above.
• A primal partial step (otherwise referable to as half step, though not necessarily having to be 50% of a full step but simply a portion of the full step) may need to be taken, however, if a primal full step would otherwise not satisfy dual feasibility. A primal no step (e.g., no primal step) is taken when the maximum number of constraints that may be active at one time is met and, thus, a constraint needs to be dropped to continue on to the next problem to be solved. For example, if thestep516 determines a primal half or full step is to be taken, the appropriate constraint is dropped instep518 to return to the same previous problem and the Lagrange multipliers and optimization variables instep514 are updated accordingly without the dropped constraint instep514. The algorithm repeats until all possible constraints are satisfied instep506 such that the algorithm proceeds to step508 to exit with the desired control outputs as a minimum feasible solution to provide to the plant instep418 of theoperation400 ofFIG. 4.
• Thus, in embodiments, and referring toFIGS. 1-5, a method for controlling theinternal combustion engine20 ofFIG. 1 that has theVGT48, theEGR valve40, and theEGR throttle44 includes solving a linear quadratic problem as describe herein with a predictive model that includes an updating algorithm to determine (i) a requested optimized engine turbine (VGT) lift66 that meets one or more active constraints, and (ii) a requested optimized EGRvalve flow rate68 that meets the one or more active constraints. Solving the linear quadratic problem includes determining, at each iteration, whether to take a determined step that is a primal partial step or a primal full step, and taking the determined step at each iteration until the linear quadratic problem is solved by the updating algorithm. Taking the primal partial step includes performing an iterative calculation, and taking the primal full step includes performing a direct calculation.
• Referring toFIGS. 1-3, the method further includes generating the requested optimizedengine turbine lift66 by controlling theVGT48 responsive to engine intake manifold pressure62 (e.g., as measured by the intake manifold pressure sensor32). For example, through control of the VGT position80 (e.g., corresponding to an angling of the turbine input vanes) by actuators, the actuators thus control theVGT48 and affect the engineintake manifold pressure62. Additionally, the method includes generating the requested optimized EGRvalve flow rate68 by controlling theEGR valve40 and theEGR throttle44 responsive to the EGR rate64 (e.g., as measured by the EGR rate sensor43). For example, through respective control of theEGR valve position82A and theEGR throttle position82B by respective actuators, the actuators thus respectively control theEGR valve40 and theEGR throttle44 and affect theEGR rate64.
EXAMPLESExample 1
• FIGS. 7A-7C show an example, along with Table 1 below, of an application of the HQPKWIK algorithm ofFIG. 5.

• TABLE 1
  		Max. Step
  	Possible	Length	Max Primal		Possible		Primal
  Iteration	Constraint	for Dual	Step	Min.	Constraint	Active	Step
  ( i)	to Add	Feasibility	Length	Length	to Drop	Set list	Taken
  0 ( 0)	4	∞	0.15	0.15	—	$\hspace{1em}\left[\begin{array}{c}0\\ 0\\ 0\\ 4\\ 0\end{array}\right]$	Full Step (Add 4)
  1 ( 1)	2	0.46	1.88	0.46	4	$\hspace{1em}\left[\begin{array}{c}0\\ 0\\ 0\\ 0\\ 0\end{array}\right]$	Half Step (Drop 4)
  2 ( 2)	2	∞	0.36	0.36	—	$\hspace{1em}\left[\begin{array}{c}0\\ 2\\ 0\\ 0\\ 0\end{array}\right]$	Full Step (Add 2)
  3 ( 3)	1	∞	0.084	0.084	—	$\hspace{1em}\left[\begin{array}{c}1\\ 2\\ 0\\ 0\\ 0\end{array}\right]$	Full Step (Add 1)
  4 ( 4)	3	0.81	∞	0.81	1	$\hspace{1em}\left[\begin{array}{c}0\\ 2\\ 0\\ 0\\ 0\end{array}\right]$	No Step (Drop 1)
  5 ( 5)	3	0.14	0.94	0.14	2	$\hspace{1em}\left[\begin{array}{c}0\\ 0\\ 0\\ 0\\ 0\end{array}\right]$	Half Step (Drop 2)
  6 ( 6)	3	∞	0.01	0.01	—	$\hspace{1em}\left[\begin{array}{c}0\\ 0\\ 3\\ 0\\ 0\end{array}\right]$	Full Step (Add 3)
  7 ( 7)	1	∞	0.003	0.003	—	$\hspace{1em}\left[\begin{array}{c}1\\ 0\\ 3\\ 0\\ 0\end{array}\right]$	Full Step (Add 1)
  8	—	—	—	—	—	—	—
  ( 8)

• The system of the example ofFIGS. 7A-7C includes two (2) optimization variables and five (5) constraints. The examples solves for the quadratic problem of Equation 3.2 and includes the following term(s), matrice(s), and vector(s):
• $Q=\left[\begin{array}{cc}\mathrm{.8814}& -\mathrm{.155}\\ -\mathrm{.155}& \mathrm{.0679}\end{array}\right]$$\mathscr{H}=\left[\begin{array}{c}\mathrm{.8109}\\ -\mathrm{.0834}\end{array}\right]$$=\left[\begin{array}{cc}-\mathrm{.7157}& -\mathrm{.6985}\\ \mathrm{.9916}& \mathrm{.1291}\\ \mathrm{.9977}& \mathrm{.0673}\\ \mathrm{.3863}& \mathrm{.9224}\\ -\mathrm{.866}& -\mathrm{.4636}\end{array}\right]$$=\left[\begin{array}{c}-3.623\\ 1.431\\ 1.202\\ 1.839\\ -4.1753\end{array}\right]$
• After the unconstrained (global) minimum

  

  ₀* is found in step502 (as shown atiteration 0 as

  

  ₀in Table 1), the example determines the most violated possible constraint instep504. In the example ofFIGS. 7A-7C, the possible constraints would be Constraints1-4 as they border a feasible region.Constraint5 would not be possible, and thus would be an inactive constraint. Further, each of the Constraints1-5 may include a single optimization variable as a constraint or multiple constraining optimization variables. For example,Constraint1 may include a temperature constraint, a pressure constraint, or a combination of a temperature or pressure constraints. Constraints may include limits such as maximum levels, minimum levels, or operating range restrictions, for example.
• As shown in the first data row of Table 1, the most violated possible constraint would beConstraint4. After determining that all possible constraints are thus not satisfied instep506, the algorithm continues on to steps510-514 to solve for a maximum step length for dual feasibility. This is found to be a maximum primal step length, as shown in the first data row, third and fourth columns, respectively, of Table 1. If the minimum length is the maximum primal step length, as found here, the algorithm continues on to take a primal step and to add the evaluated constraint (i.e., Constraint4) to an active set list, as shown in Table 1, and thus followssteps516 and520 to return tostep504.
• Atiteration 1, where the system is now at point

  

  ₁along a boundary represented byConstraint4, the algorithm determines instep504 that the next most violated possible constraint is now Constraint2 (as shown in the second data row of Table 1). As steps510-514 are applied toConstraint2, however, the maximum step length for dual feasibility is found to be less than the maximum primal step length. This indicates that the full primal step would be out of an error boundary zone and infeasible with bothConstraints2 and4, and thusConstraint4 needs to be dropped. A half (e.g., partial) step is then taken atstep516 to point

  

  ₂,Constraint4 is dropped instep518, and the Lagrange multipliers and optimization variables are updated atstep514 for the previous same problem withoutConstraint4. In iteration 2 (with the algorithm now starting at point

  

  ₂),Constraint2 is again evaluated (this time without Constraint4) and the Lagrange multipliers and optimization variables are updated atstep514. Atstep516, the maximum primal step length and the maximum step length for dual feasibility are compared, and the maximum primal step length is found to be the minimum length. Thus,step516 determines that a full (primal) step should be taken and thatConstraint2 should be added to the active set list to arrive at point

  

  ₃. The algorithm continues throughiterations 3 through 8 until all possible constraints are satisfied atstep506.
• More particularly, atiteration 4 when the algorithm is at point

  

  ₄and attempts to addConstraint3, a determination is made atstep516 that no step should be taken as a maximum number of constraints are active at one time. Thus, the algorithm updates to point

  

  ₅that is the same as point

  

  ₄, andConstraint1 is dropped. Initeration 5, however, it is determined atstep516 to take a half step to point

  

  ₆, andConstraint2 is dropped from the active set list. At iteration 6, it is determined atstep516 that a full step may be taken to addConstraint3 instep520. In iteration 7,Constraint1 is added in a full step moving the algorithm from point

  

  ₇, to point

  

  ₈, which is a most minimum point of the feasible region that satisfies all possible restraints and is an output solution for the plant (e.g., utilized instep418 ofFIG. 4).
Example 2
• Referring toFIGS. 8A-8B, another example is presented with respect to use of the algorithm ofFIG. 5.FIG. 8A shows a graphical representation of the primal problem having two (2) optimization variables and four (4) constraints. For example, for a given

  

  and

  

  as set forth below, the problem involves an optimal value of a VGT position

  

  ₁* (units of % closed) and an EGR valve opening

  

  ₂* (units of % open) that minimizes the quadratic cost function J(

  

  ) and satisfies the following four constraints of Table 2:

• TABLE 2
  		Mathematical
  ConstraintNumber	Description	Representation
  1	Minimum VGT % Closed =	h1(   ) = − 1+ 40 ≤ 0
  	40%
  2	Maximum VGT %	h2(   ) =  1− 80 ≤ 0
  	Closed = 80%
  3	Minimum EGR Valve	h3(   ) = − 2+ 2 ≤ 0
  	Opening = 2%
  4	Maximum EGR Valve	h4(   ) =  2− 70 ≤ 0
  	Opening = 70%

• In particular, the associated primal problem is:
• $\underset{}{\min }J\left(\right)=\frac{1}{2}{}^T+{\mathscr{H}}^T$$s.t.h\left(\right)=\mathrm{}+\le 0$$\mathrm{where}=\left[\begin{array}{c}{}_1\\ {}_2\end{array}\right],Q=\left[\begin{array}{cc}10.7& 0\\ 0& 26\end{array}\right],\mathscr{H}=\left[\begin{array}{c}-321\\ -26\end{array}\right],=\left[\begin{array}{cc}-1& 0\\ 1& 0\\ 0& -1\\ 0& 1\end{array}\right],\mathrm{and}=\left[\begin{array}{c}40\\ -80\\ 2\\ -70\end{array}\right].$
• TheHQPKWIK algorithm500 first starts atstep501 and then sets the optimization variables to an unconstrained minimum atstep502 such that, as shown inFIG. 8B,
• ${}_0=-Q^{-1}\mathscr{H}=\left[\begin{array}{c}30\\ 1\end{array}\right].$
• Atstep504, the most violated constraint is determined to beConstraint1, as max(

  

  −

  

  >0)=Constraint1. Instep506, thealgorithm500 determines that not all the possible constraints are satisfied (as it had identifiedConstraint1 as a most violated possible constraint in step504), and thealgorithm500 proceeds to step510, in which the active set list is converted to a unique identifier. In the present example, the active set list is of the form S={Constraint4,Constraint3,Constraint2, Constraint1} such that S₁={0,0,0,1}⇒0001_bin=1_dec. Thealgorithm500 then proceeds to step512 to calculate primal and dual step direction and determine step length, which it does to find that a full primal step should be taken andConstraint1 added as t=min(10,∞)=10 (primal space)⇒Full Primal Step. Thealgorithm500 proceeds to step514 to update the optimization variables

  

  ₁and the Lagrange multipliers μ₁as follows:
• ${}_1=-{\mathrm{<figure-callout\; label="z"\; filenames="US20190107072A1-20190411-D00006.png,US20190107072A1-20190411-D00007.png"\; state="\{\{state\}\}">z</figure-callout>}}_{1\mathrm{dec}}\mathscr{H}+{\mathrm{<figure-callout\; label="r"\; filenames="US20190107072A1-20190411-D00006.png,US20190107072A1-20190411-D00007.png"\; state="\{\{state\}\}">r</figure-callout>}}_{1\mathrm{dec}}^T{}_1=\left[\begin{array}{c}40\\ 2\end{array}\right]$${\mathrm{<figure-callout\; label="\mu "\; filenames="US20190107072A1-20190411-D00006.png,US20190107072A1-20190411-D00007.png"\; state="\{\{state\}\}">\mu </figure-callout>}}_1=0-\left[\begin{array}{c}-{\mathrm{<figure-callout\; label="r"\; filenames="US20190107072A1-20190411-D00006.png,US20190107072A1-20190411-D00007.png"\; state="\{\{state\}\}">r</figure-callout>}}_{1\mathrm{<figure-callout\; label="dec\; \; v"\; filenames="US20190107072A1-20190411-D00006.png,US20190107072A1-20190411-D00007.png"\; state="\{\{state\}\}">dec</figure-callout>}}v_{}{}_1\\ 1\end{array}\right]t=\left[\begin{array}{c}10\\ 1\end{array}\right].$
• Afterwards, thealgorithm500 proceeds to step516 as a full primal step is taken (between points

  

  ₀and

  

  ₁as shown inFIG. 8C) andConstraint1 is added before thealgorithm500 returns to step504 to re-start the process with the newly addedConstraint1.
• For example, instep504, thealgorithm500 determines the new most violated (possible) constraint to beConstraint3, as max(

  

  −

  

  >0)=Constraint3 and thus notes instep506 that not all (possible) constraints are satisfied, moving on to step510. Instep510, the active set list is converted to a unique identifier as follows: S_1,3={0,1,0,1}⇒0101_bin=5_dec. The dual and primal space arrays associated with theunique identifier5 is retrieved, and, instep512, the primal and dual step directions are calculated to determine step length to find that t=min(1,∞)=1 (primal space)⇒Full Primal Step. Thealgorithm500 proceeds to step514 to update the optimization variables

  

  ₂and the Lagrange multipliers μ₂as follows:
• ${}_2=-{\mathrm{<figure-callout\; label="z"\; filenames="US20190107072A1-20190411-D00006.png,US20190107072A1-20190411-D00007.png"\; state="\{\{state\}\}">z</figure-callout>}}_{5\mathrm{dec}}\mathscr{H}+{\mathrm{<figure-callout\; label="r"\; filenames="US20190107072A1-20190411-D00006.png,US20190107072A1-20190411-D00007.png"\; state="\{\{state\}\}">r</figure-callout>}}_{5\mathrm{dec}}^T{}_2=\left[\begin{array}{c}40\\ 2\end{array}\right]$${\mathrm{<figure-callout\; label="\mu "\; filenames="US20190107072A1-20190411-D00005.png,US20190107072A1-20190411-D00006.png"\; state="\{\{state\}\}">\mu </figure-callout>}}_2={\mathrm{<figure-callout\; label="\mu "\; filenames="US20190107072A1-20190411-D00006.png,US20190107072A1-20190411-D00007.png"\; state="\{\{state\}\}">\mu </figure-callout>}}_1+\left[\begin{array}{c}-{\mathrm{<figure-callout\; label="r"\; filenames="US20190107072A1-20190411-D00006.png,US20190107072A1-20190411-D00007.png"\; state="\{\{state\}\}">r</figure-callout>}}_{5\mathrm{<figure-callout\; label="dec\; \; v"\; filenames="US20190107072A1-20190411-D00005.png,US20190107072A1-20190411-D00006.png"\; state="\{\{state\}\}">dec</figure-callout>}}v_{}{}_2\\ 1\end{array}\right]t=\left[\begin{array}{c}10\\ 1\end{array}\right].$
• Afterwards, thealgorithm500 proceeds to step516 as a full primal step is taken (between points

  

  ₁and

  

  ₂as shown inFIG. 8D) andConstraint3 is added before thealgorithm500 returns to step504 to re-start the process with the newly addedConstraint3. This time, atstep504, the process determines that there is not another most violated possible constraint (as max(

  

  −

  

  >0)=Null), and then, instep506, all possible constraints are determined to be satisfied at point

  

  ₁(such that thealgorithm500 exits atstep508 with a solution of

  

  ₂). As shown inFIG. 8D, point

  

  ₂is the minimum point or solution at the edge of the primal feasible region that is constrained by the four constraints.
• Referring toFIGS. 9-10, to validate the HQPKWIK algorithm, drive cycle data for a diesel engine with a linear MPC controlling the air path (e.g., the plant) was collected as transient bench data for a WLTP cycle and was evaluated. Additionally, the HQPKWIK algorithm as described herein was used as the solver and evaluated. To validate the results, each individual problem corresponding to each time step was given to a Matlab solver (offline), and a maximum error between the HQPKWIK algorithm and transient bench data results was determined, as shown inFIG. 10.
• Accumulation errors due to, for example, an iterative approach may be larger than an acceptable value of about 1e-6, indicating the errors are larger than a bound and the algorithm may not work as expected. By utilizing the HQPKWIK algorithm and a hybrid approach of utilizing iterative and direct calculations in primal space, the errors may be maintained below the acceptable value such that evaluated constraints are satisfied, as shown inFIG. 10. The results ofFIG. 10 were within an indicative error range of 1.5e-5%, illustrating a successful ability to calculate the solution via the HQPKWIK algorithm correctly online.
• The HQPKWIK algorithm is a hybrid range space algorithm that allows for a constrained optimization solver for production use in a plant of a diesel engine air path to be applied to a problem that changes with time (time-varying) while using an updating algorithm. The updating algorithm updates optimization variables and Lagrange multipliers of the problem based on a step length that is determined to be taken in primal and dual spaces. The HQPKWIK algorithm also successfully meets challenges such as working within the confines of a low ROM size (i.e., 4.63 KB total), being able to run on an ECU, and being able to find a correct solution accurately and a finite and quick timeframe. For example, as set forth above,FIG. 10 shows a maximum error of about 1.5e-5%, and the run example resulted in an average of about 250 FLOPS. To note, the computation time and ROM size grows as more constraints are added as HQPKWIK is an active set method, and thus the solver described herein is generally useful for solving any constrained optimization convex problem with a small number of constraints that have linear inequality.
• It is noted that the terms “substantially” and “about” and “approximately” may be utilized herein to represent the inherent degree of uncertainty that may be attributed to any quantitative comparison, value, measurement, or other representation. These terms are also utilized herein to represent the degree by which a quantitative representation may vary from a stated reference without resulting in a change in the basic function of the subject matter at issue.
• While particular embodiments have been illustrated and described herein, it should be understood that various other changes and modifications may be made without departing from the spirit and scope of the claimed subject matter. Moreover, although various aspects of the claimed subject matter have been described herein, such aspects need not be utilized in combination. It is therefore intended that the appended claims cover all such changes and modifications that are within the scope of the claimed subject matter.

Claims (10)

What is claimed is:

1. A method for controlling an internal combustion engine having a variable geometry turbine (VGT), an exhaust gas recirculation (EGR) valve, and an EGR throttle, the method comprising:

solving a linear quadratic problem with a predictive model comprising an updating algorithm in order to determine: (i) a requested optimized VGT lift that meets one or more constraints; and (ii) a requested optimized EGR valve flow rate that meets the one or more constraints, wherein solving the linear quadratic problem includes:

determining whether to take a determined step comprising one of a primal partial step and a primal full step at each iteration, and

taking the determined step at each iteration until the linear quadratic problem is solved by the updating algorithm, wherein:

taking the primal partial step comprises performing an iterative calculation, and

taking the primal full step comprises performing a direct calculation;

generating the requested optimized VGT lift responsive to an engine intake manifold pressure by controlling the VGT; and

generating the requested optimized EGR valve flow rate responsive to an EGR rate by controlling the EGR valve and the EGR throttle.

2. The method ofclaim 1, further comprising:

setting one or more engine operating parameters as the one or more constraints to form a non-linear problem; and

deriving the linear quadratic problem based on the non-linear problem, wherein the linear quadratic problem is convex and time-varying.

3. The method ofclaim 2, wherein the linear quadratic problem to solve comprises the following equation:

$\underset{}{\min }J\left(\right)=\frac{1}{2}{}^T\mathrm{}+{\mathscr{H}}^T$$s.t.\mathrm{}\le $

where

is a time varying linear term in the cost function;

is a constant quadratic term in the cost function;

is a constant constraint matrix;

is a time varying constraint vector; and

is one or more optimization variables comprising time varying control inputs.

4. The method ofclaim 3, wherein the time varying control inputs comprise at least one of the a VGT lift and an EGR flow, and further comprising adding Lagrangian multipliers μ to determine a solution for optimization variable(s)

* and Lagrangian multiplier(s) μ* that are greater than or equal to zero to satisfy dual feasibility, while utilizing only active, feasible constraints of the one or more constraints.

5. The method ofclaim 4, wherein:

the VGT comprises turbine input vanes configured to be angled to be opened, partially opened, or closed; and

the active, feasible constraints comprises at least one of a maximum EGR rate, a minimum EGR rate, a maximum EGR flow command, a minimum EGR flow command, a maximum VGT lift closed command to control a maximum amount of closure of the turbine input vanes, and a minimum VGT lift closed command to control a minimum amount of closure of the turbine input vanes.

6. The method ofclaim 5, further comprising simplifying the determined solution for a current active set list L and to define iterative and direct calculation approaches to determine how to move from the current active set list L to a new active set list L+1, resulting in the following simplified HQPKWIK equation:

${}_{L+1}=\{\begin{array}{cc}{}_L+{\mathrm{zv}}_{\mathrm{knext}}t& \mathrm{Primal}\mathrm{Partial}\left(\mathrm{Half}\right)\mathrm{Step}\\ -{T_{L+1}^U\left(T_{L+1}^U\right)}^T\mathscr{H}+T_{L+1}^CR_{L+1}^{-T}{}_{L+1}& \mathrm{Primal}\mathrm{Full}\mathrm{Step}\end{array}\begin{array}{cc}{\mu }_{L+1}={\mu }_L-{\mathrm{rv}}_{\mathrm{knext}}t& \mathrm{Dual}\mathrm{Step}\end{array}$

where t is a minimum length that still maintain dual feasibility, v_knextis a next constraint to be added to a current constraint set list, z is an array that is indicative of a search direction in primal space, and r is an array that is indicative of a search direction in dual space, such that:

z=T_L^U(T_L^U)^T

r=R_L^{7 −1}(T_L^C)^T.

7. The method ofclaim 6, further comprising computing the z and r arrays on-line.

8. The method ofclaim 6, further comprising pre-computing portions of the z and r arrays off-line.

9. The method ofclaim 8, further comprising retrieving the z and r arrays after a determination of a unique identifier at least partially based on the current active set list.

10. The method ofclaim 9, wherein the determination of the unique identifier utilizes a conversion of one or more binary numbers of the current active set list to a decimal value, and the unique identifier comprises the decimal value.

Images