Control of thermal conditions in air-conditioned environments based on multivariate optimization

Abstract

To provide a system and a method for controlling an operation of HVAC system.SOLUTION: A system includes at least one processor and a memory. The memory stores instructions that, when executed by the at least one processor: causes the system to receive data 100 indicative of an operation of HVAC system and optimize a cost function to reduce a difference between a full thermal state 104 of a conditioned environment simulated from the received data on the basis of a physical model of airflow including a partial differential equation (PDE) and a thermal state 112 of the conditioned environment reconstructed by an observer using a reduced order model (ROM) of airflow including an ordinary differential equation (ODE) from values of the full thermal state at a set of locations. Therein, optimization of the cost function is a multivariable optimization of a structure of the observer and a structure of the ODE.SELECTED DRAWING: Figure 1

Description

[Technical field]

[0001] FIELD OF THE DISCLOSURE This disclosure relates generally to heating, ventilation and air conditioning (HVAC) systems, and more particularly to controlling the operation of a heating, ventilation and air conditioning (HVAC) system arranged to regulate an environment. [Background technology]

[0002] background Air conditioning systems, such as heating, ventilation and air conditioning (HVAC) systems, are widely used in industrial and domestic applications to regulate the environment. The introduction of variable speed compressors, variable position valves and variable speed fans in the vapor compression cycle has greatly increased the operational flexibility of HVAC systems. The efficiency of an HVAC system can be improved by controlling the variable components of the HVAC system. In other words, control inputs can be given to the HVAC system to improve its efficiency. For example, the compressor speed can be adjusted to regulate the refrigerant flow rate. The evaporator and condenser fan speeds can be changed to vary the heat transfer between the air and the heat exchanger. Changing the opening of the expansion valve can affect the pressure drop between high and low pressure in the HVAC system, affecting the refrigerant flow rate and the corresponding superheat temperature at the evaporator outlet.

[0003] The combination of control inputs to an HVAC system that delivers a given amount of heat is often not unique, and various combinations of control inputs consume different amounts of energy. Therefore, it is desirable to operate the HVAC system using a combination of control inputs that minimizes energy and maximizes the energy efficiency of the HVAC system. Some methods for maximizing energy efficiency rely on the use of mathematical models of the physics of the HVAC system. Such model-based methods attempt to describe the effects of HVAC system control inputs on the thermodynamic behavior and energy consumption of the HVAC system. In such model-based methods, a pre-analytical model is used to predict the combination of control inputs that will meet the heat load requirements and minimize energy.

[0004] However, the use of mathematical models to predict such control input combinations has drawbacks. For example, the mathematical models rely on simplifying assumptions to generate mathematically tractable representations. Such assumptions ignore important effects or do not consider installation-specific characteristics such as room size, causing the mathematical models to deviate from the actual behavior of the HVAC system. In some schemes, knowledge of the airflow that occurs during operation of the HVAC system is used to optimize the operation of the HVAC system by formulating a physical model of the airflow. However, the physical model of the airflow is infinite-dimensional and very complex to use in real-time control applications. Also, the physical model of the airflow may change during the operation of the air conditioning system.

[0005] Prediction and control of the Boussinesq equations (which consist of the Navier-Stokes equations together with the energy conservation equation), which are a physical model of airflow, are computationally expensive, making their application to real-time control difficult. It is therefore desirable to rapidly reconstruct the humidity, temperature and velocity of airflow in an environment from sensory data. To accomplish this task, various observer design techniques can be employed, such as Kalman filters, H-∞ filters, etc. Such observer design techniques are based on reducing the estimation error in a statistical sense for a given norm that measures the estimation error. The estimation error can be defined as the norm of the difference between the reconstructed quantities and the ground truth (e.g., a CFD solution). Once an observer is designed, the quantities of interest (e.g., air humidity, temperature and velocity) in all or part of the environment (e.g., an occupied area) can be reconstructed and then used for control (e.g., closed-loop control) to achieve a desired state (e.g., thermal comfort of occupants or minimum energy consumption of an HVAC system).

[0006] However, reconstruction / estimation of the quantities of interest and real-time control of airflow are challenging. For example, due to the spatiotemporal evolution of the quantities of interest, the location of the sensor plays a key role in the estimation. If the sensor is located far from the inlet of the HVAC system, the estimation of the quantities of interest near the inlet will be poor or impossible. Furthermore, the environment often contains multiple coherent structures. If the sensor is placed far from such structures, for example in the region of a large air vortex that occurs between the inlet and the exhaust of the HVAC system, the reconstruction error will be large. The number and location of sensors required to accurately and robustly reconstruct the quantities of interest may vary depending on the relative positions of the inlets and outlets, the number of inlets and outlets (actuators), heat losses through walls, and the configuration of objects and occupants in the environment, among other factors.

[0007] Furthermore, there is a great deal of research into selecting scientifically correct and experimentally verified models of airflow in the environment. However, such models are partial differential equations (PDEs), for which analytical solutions often do not exist. Thus, numerical solutions are required. Such numerical solutions (e.g., finite difference, finite volume, finite element) are tedious and require huge computational costs and hardware resources. This is prohibitive for real-time control applications. Therefore, the use of reduced order models (ROMs) is desirable. The advantage of using ROMs rather than the original governing equations is that they have a smaller order than the discretized algebraic equations resulting from PDEs, making them applicable to real-time control applications. However, various ROMs cannot inherit the original system characteristics, for example, in terms of the stability of HVAC systems. Summary of the Invention

[0008] overview An objective of some embodiments is to determine a stable reduced order model (ROM) that can be used to estimate airflow dynamics during operation of a heating, ventilation, and air conditioning (HVAC) system. Additionally, an objective of some embodiments is to use the determined ROM for real-time control of operation of the HVAC system such that the thermal load requirements of the operation are met and the performance of the HVAC system is optimized. Additionally or alternatively, an objective of some embodiments is to determine sensor locations along with the ROM.

[0009] The HVAC system may be arranged to condition an environment. The environment may be a room or space of a building, or the entire building in which the HVAC system is installed. In some embodiments, the environment may correspond to a space of a building in which an occupant is located or occupied. The HVAC system is configured to output air in the environment to condition the environment and ensure thermal comfort of the occupant of the environment. The thermal condition of the output air includes the temperature and velocity of the air output by the HVAC system. In an alternative embodiment, the thermal condition of the output air includes one or a combination of the temperature, velocity, and humidity of the air output by the HVAC system.

[0010] According to an embodiment, a system for controlling operation of an HVAC system is provided. The system may include a processor and a memory configured to store instructions executable by the processor. The processor is configured to execute the instructions stored in the memory to perform steps necessary to control operation of the HVAC system. The system may control the HVAC system to optimize energy consumption of the HVAC system by varying thermal conditions in an environment. In an embodiment, the system is configured to receive data indicative of operation of the HVAC system. The data indicative of operation of the HVAC system may include boundary conditions corresponding to the environment. Additionally, the data indicative of operation of the HVAC system may include a geometry of the environment, a number of inlets and outlets, etc.

[0011] Since HVAC control optimizes energy consumption by changing the thermal state in the environment, it is desirable to know the dynamics of the airflow that occurs in the environment during the operation of the HVAC system. For this purpose, a physical model of the airflow is considered. For example, the physical model of the airflow can be defined according to the Boussinesq equation. However, the computational cost of predicting and controlling the Boussinesq equation (which consists of the Navier-Stokes equation together with the energy conservation equation), which is the physical model of the airflow, is difficult for real-time control (or online control) of the HVAC system. For example, the Boussinesq equation, which describes the physical model of the airflow, is a partial differential equation (PDE), and it is difficult to solve it in real time. Therefore, it is desirable to express the physical model of the airflow including the PDE as a reduced-order model (ROM) represented by an ordinary differential equation (ODE) with fewer parameters than the number of parameters of the physical model of the airflow.

[0012] Some embodiments are based on the understanding that many different techniques can perform such dimensionality reduction. For example, based on data describing the operation of an HVAC system and the Boussinesq equation, the system performs a computational fluid dynamics (CFD) simulation to generate an estimate of the thermal state at every point (or location) of the environment. That is, the CFD simulation technique can be used to simulate the airflow dynamics in the environment. The thermal state at every point generated by the CFD simulation is called the full thermal state. The full thermal state is considered as the ground truth value of the thermal state at every point of the environment. Furthermore, a ROM is derived such that using the ROM to estimate the airflow dynamics in the environment approximates the full thermal state.

[0013] For example, in a CFD simulation, a full thermal state occurs at each instant in time. The full thermal state at an instant in time is called a snapshot. The different snapshots of the CFD simulation can be organized into a matrix that captures the airflow dynamics, and the ROM is determined by solving an optimization problem such that the difference between the airflow dynamics estimated by the ROM and the airflow dynamics provided by the CFD simulation is minimized in a low-dimensional space. Such a scheme is called the minimum ||AX i -X i+1 It can be thought of as a linear estimation of the ROM as a matrix A by solving the minimization problem given by ||, where matrix A is the ROM and X i represents the thermal state at time i. The matrix A is used to linearly estimate the airflow dynamics in a low-dimensional space (e.g., x i+1 =A * x i ), can be used. This minimization problem can be solved by several matrix decomposition techniques. For example, matrix A can be determined using Proper Orthogonal Decomposing (POD) with Galerkin Projection, Dynamic Mode Decomposition (DMD), Spectral POD, and Balanced POD. Using the techniques described above, the minimization problem can be solved and matrix A, i.e., ROM, can be determined.

[0014] In some cases, the determined ROM may not be linear. The determined ROM may be a quadratic equation or a nonlinear ROM. Moreover, the above-mentioned techniques capable of solving the minimization problem, alone or in combination, do not guarantee that the determined ROM can be used for online HVAC control, for example, the determined ROM may lack stability or observability. Since the ROM operates in a low-dimensional space, an observer is needed that can generate an estimate / reconstruction of the thermal state at all desired points in the environment using the airflow dynamics captured by the ROM in the low-dimensional space and actual measurements at several points in the environment. According to an embodiment, the observer is a function of the ROM, i.e., the ODE and / or matrix A, and the placement of the sensors in the environment. Specifically, the structure of the observer includes a gain, referred to herein as the gain matrix K, and designing the observer means determining the matrix K.

[0015] Some embodiments are based on the recognition that since matrix K depends on matrix A, matrix A is determined first, and then matrix K is determined based on matrix A. However, some embodiments are based on the recognition that the calculation of matrix K may affect the determination of matrix A. In other words, there is an interdependency between the determination of matrix A and the determination of matrix K. As a result, some embodiments are based on the recognition that the joint and interdependent optimization of matrix A and matrix K may improve the accuracy and global optimality of the reconstruction of the thermal state. To that end, the system performs a multivariate optimization including the joint optimization of the structure of the observer (i.e., matrix K) and the structure of the ROM (i.e., matrix A). That is, the ROM and the observer are determined together as optimization parameters of a cost function that uses the complete thermal state as the ground truth.

[0016] Therefore, the system optimizes a cost function to reduce the difference between the complete thermal state of the environment simulated based on the received data and a physical model of the airflow including the PDE and the thermal state of the environment reconstructed by the observer using the ROM from the values ​​of the complete thermal state at a set of positions. The optimization of the cost function is a multivariate optimization of the structure of the observer and the structure of the ODE or ROM. According to an embodiment, a set of positions of sensors for measuring the thermal state is predetermined.

[0017] Additionally or alternatively, some embodiments are based on the recognition that matrix K is a function of the placement of sensors in the environment. The sensors may measure one or a combination of temperature, velocity, and humidity of air in the environment. The placement of sensors can be represented as a binary matrix C with elements corresponding to locations in the environment. The elements of matrix C have a value of 1 in the elements corresponding to the locations of the sensors and a value of 0 otherwise. In some embodiments, matrix C, i.e., the locations of the sensors, are predetermined by the designer of the HVAC system. However, in some other embodiments, matrix C may be unknown and / or may be a design parameter. Some embodiments are based on the recognition that matrix K depends on matrices A and C, and thus matrix C may be incorporated as an optimization parameter in a multivariable optimization. In other words, the set of ODE or ROM structure, observer structure, and sensor locations are jointly optimized as part of a multivariable optimization. This optimization is called "co-design." Co-design not only results in an accurate reconstruction of the thermal state, but also optimizes the sensor locations.

[0018] Some embodiments are based on another realization that when sensor placement is also part of the multivariable optimization, it is computationally efficient to perform the multivariable optimization via pairwise combinations such that an observer is determined considering an ODE and a set of sensor locations, an ODE is determined considering an observer and a set of sensor locations, and a set of sensor locations is determined considering an observer and an ODE. In this way, the convergence of the multivariable optimization can be improved.

[0019] Furthermore, according to some embodiments, the controller of the HVAC system can be modified with the optimized structure of the observer and the optimized structure of the ODE.

[0020] According to some embodiments, in multivariate optimization, the optimization of both the observer structure and the ODE structure are the primary optimization objectives. In other words, in multivariate optimization, the observer and the ODE are given equal / balanced importance. However, some embodiments are based on the recognition that the importance of the ODE and the observer in the objective of online HVAC control is not equal. In some embodiments, since the objective of the observer is to enable the use of the ODE, the ODE is the primary objective and the observer is a secondary objective, i.e., the observer is subordinate to the ODE. Therefore, only considering the ODE and the observer jointly may not be optimal.

[0021] To alleviate such problems, the multivariable optimization is transformed into a multivariable constrained optimization, which achieves an interdependent yet unbalanced multivariable optimization of the ODE and the observer, with the ODE as the first optimization objective and the observer as the second optimization objective. That is, in the multivariable constrained optimization, the ODE is weighted more than the observer. Thus, the ODE and the observer are given unequal / disproportionate importance. To achieve this interdependence, the system optimizes the structure of the observer and the structure of the ODE jointly, e.g., as part of one unit multivariable optimization. To achieve the disproportionate importance, in some embodiments, the multivariable constrained optimization is used in preference to the optimization of the ODE, where the optimization of the observer structure acts as a constraint on the optimization of the ODE. In this way, the optimization of the observer structure becomes subordinate to the optimization of the ODE, since the observer structure does not have an independent optimization objective, but serves only to constrain the structure of the ODE.

[0022] Furthermore, some embodiments are based on the recognition that when multivariate optimization is transformed into multivariate constrained optimization with unbalanced optimization objectives, multivariate constrained optimization can be used to address some potential deficiencies in the structure of ODE and observer. For example, some approaches that sequentially determine ODE and observer, i.e., first determine ODE and then determine observer considering ODE, may make matrix A of ODE unstable. As used herein, stability is the property that a solution that starts near an equilibrium point will stay near the equilibrium point forever. Therefore, even if the PDE is unstable, it is still necessary to determine a stable ROM for control / observation.

[0023] Some embodiments are based on the recognition that because multivariate optimization adds flexibility to the optimization of the structure of the ODE, it may be advantageous to perform a multivariate optimization subject to a stability constraint imposed on the structure of the ODE. Such optimization is sometimes referred to as constrained optimization. According to one embodiment, the stability constraint may be a Lyapunov inequality.

[0024] Additionally or alternatively, some embodiments are based on the recognition that constraints imposed in constrained optimization may include constraints on the structure of the observer without significantly increasing the computational complexity of the multivariate optimization. Examples of such constraints include observability constraints. As referred to herein, observability is a measure of how well the internal state of a system, e.g., the temperature of the environment, can be inferred from knowledge of its external outputs, e.g., the measurements of a thermostat. Enforcing observability on the structure of the observer may therefore be advantageous to ensure that the temperature in the environment can be estimated with a given set of measurements. Examples of observability constraints include the Lyapunov equation of the observability Gramian, with the observability matrix as the maximum rank. Such constrained optimization results in a stable observable ROM.

[0025] However, some embodiments are based on the recognition that the conversion of ODE optimization to multivariate optimization, and in particular to multivariate constrained optimization, introduces additional problems that are not present in the usual single-variable optimization of the structure of ODEs. For example, the aforementioned techniques that can solve the optimization of the structure of ODEs are not possible or practical for multivariate constrained optimization. According to some embodiments, the multivariate constrained optimization can be solved using semidefinite programming (SDP). SDP is a subdiscipline of convex optimization that concerns the optimization of linear objective functions on the intersection of a cone of semidefinite matrices with an affine space. The system optimizes the cost function using SDP, which estimates the structure of an observer and the structure of an ODE that are dependent on each other.

[0026] Furthermore, some embodiments are based on the recognition that SDP is a sub-discipline of convex optimization, and although multivariate optimization of observer-ODE structures is a convex problem, the dimensionality of the ODE still imposes a potentially infeasible computational burden on SDP in convex domains. However, some embodiments are based on the recognition that large convex but constrained problems are harder to solve than large non-convex but unconstrained problems, because SDP may not scale well with the large dimensions of matrix A (even though the dimensions of matrix A are significantly lower than those of CFD). To that end, some embodiments disclose an indirect SDP method for solving large unconstrained non-convex problems by converting constrained convex multivariate optimization to unconstrained non-convex multivariate optimization using Cholesky decomposition.

[0027] Accordingly, an embodiment discloses a system for controlling operation of an HVAC system arranged to condition an environment, the system comprising at least one processor and a memory, the memory storing instructions that, when executed by the at least one processor, cause the system to receive data indicative of operation of the HVAC system, and optimize a cost function from the received data to reduce a difference between a full thermal state of the conditioned environment simulated based on a physical model of airflow including partial differential equations (PDEs) and a thermal state of the conditioned environment reconstructed by an observer from values ​​of the full thermal state at a set of locations using a reduced order model of airflow (ROM) including ordinary differential equations (ODEs), the optimization of the cost function being a multivariate optimization of a structure of the observer and a structure of the ODE, and further causing the system to modify a controller of the HVAC system with the optimized structure of the observer and the optimized structure of the ODE.

[0028] Accordingly, another embodiment discloses a method for controlling operation of an HVAC system arranged to condition an environment, the method comprising: receiving data indicative of operation of the HVAC system; optimizing a cost function from the received data to reduce a difference between a thermal state of the conditioned environment simulated based on a physical model of airflow including partial differential equations (PDEs) and a full thermal state of the conditioned environment reconstructed by an observer from values ​​of the full thermal state at a set of locations using a reduced order model of airflow (ROM) including ordinary differential equations (ODEs), the optimization of the cost function being a multivariate optimization of a structure of the observer and a structure of the ODE, the method further comprising modifying a controller of the HVAC system with the optimized structure of the observer and the optimized structure of the ODE. [Brief description of the drawings]

[0029] [Figure 1]1 is a schematic diagram illustrating principles used by some embodiments to control the operation of a heating, ventilation, and air conditioning (HVAC) system arranged to condition an environment. [Figure 2A] FIG. 2 is a schematic diagram illustrating the conversion of a multivariate optimization of an observer structure and an ordinary differential equation (ODE) structure to a multivariate constrained optimization according to some embodiments. [Figure 2B] FIG. 1 is a schematic diagram illustrating constrained optimization for determining a stable reduced-order model (ROM), according to some embodiments. [Figure 2C] FIG. 2 is a schematic diagram illustrating constrained optimization subject to stability and observability constraints according to some embodiments. [Diagram 3] FIG. 1 is a schematic diagram for solving multivariate constrained optimization based on semidefinite programming (SDP), according to some embodiments. [Figure 4] FIG. 1 is a block diagram illustrating a system for controlling the operation of an HVAC system arranged to condition an environment, according to some embodiments. [Diagram 5] FIG. 1 is a block diagram for computational fluid dynamics (CFD) or partial differential equation (PDE) numerical solution according to some embodiments. [Figure 6] FIG. 1 is a schematic diagram of a co-design for determining a stable observable ROM with sensor positions according to some embodiments. [Figure 7] FIG. 2 is a block diagram illustrating a memory-limited BFGS (L-BFGS or LM-BFGS) scheme according to some embodiments. [Figure 8] 1 is a flowchart illustrating a method for controlling operation of an HVAC system arranged to condition an environment, according to some embodiments. [Figure 9] FIG. 1 illustrates control of an HVAC system using the system according to some embodiments. DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0030] Detailed Description In the following description, for purposes of explanation, numerous specific details are set forth in order to provide a thorough understanding of the present disclosure. However, it will be apparent to those skilled in the art that the present disclosure may be practiced without these specific details. In other instances, the devices and methods are shown only in block diagram form in order not to obscure the present disclosure.

[0031] As used herein and in the claims, the terms "for example" and "such as," as well as the verbs "comprise," "have," "include," and other verb forms thereof, when used in connection with a list of one or more components or other items, are to be construed as open-ended, meaning that the list is not to be considered as excluding other additional components or items. The term "based on" means based at least in part on. Furthermore, it is to be understood that the phraseology and terminology used herein are for purposes of this specification and should not be considered limiting. Any headings used herein are for convenience only and have no legal or limiting effect.

[0032] FIG. 1 is a schematic diagram illustrating principles used by some embodiments to control the operation of a heating, ventilation, and air conditioning (HVAC) system arranged to condition an environment. The environment may be a room or space of a building, or an entire building in which the HVAC system is installed. In some embodiments, the environment may correspond to a space of a building in which an occupant is located or occupied. In various implementations, the HVAC system may include multiple HVAC units installed in the environment. The HVAC system is configured to output air in the environment to condition the environment and ensure thermal comfort of the occupants of the environment. The thermal condition of the output air includes the temperature and velocity of the air output by the HVAC system. In an alternative embodiment, the thermal condition of the output air includes one or a combination of the temperature, velocity, and humidity of the air output by the HVAC system.

[0033] According to an embodiment, a system for controlling operation of an HVAC system is provided. The system may include a processor and a memory configured to store instructions executable by the processor. The processor is configured to execute the instructions stored in the memory to perform steps necessary to control operation of the HVAC system. The system may control the HVAC system to optimize energy consumption of the HVAC system by varying thermal conditions in an environment. In an embodiment, the system is configured to receive data 100 indicative of operation of the HVAC system. The data 100 indicative of operation of the HVAC system may include boundary conditions corresponding to the environment. Additionally, the data 100 indicative of operation of the HVAC system may include the geometry of the environment, the number of inlets and outlets, etc.

[0034] Since HVAC control optimizes energy consumption by changing the thermal state in the environment, it is desirable to know the dynamics of airflow occurring in the environment during the operation of the HVAC system. To do so, a physical model of airflow is considered. For example, the physical model of airflow can be defined according to the Boussinesq equation. However, the computational cost for predicting and controlling the Boussinesq equation (which consists of the Navier-Stokes equation together with the energy conservation equation), which is the physical model of airflow, is difficult for real-time control (or online control) of the HVAC system. For example, the Boussinesq equation, which describes the physical model of airflow, is a partial differential equation (PDE), and it is difficult to solve it in real time. Therefore, it is desirable to express the physical model of airflow including PDE as a reduced order model (ROM) represented by an ordinary differential equation (ODE) with fewer parameters than the number of parameters of the physical model of airflow.

[0035] Some embodiments are based on the understanding that many different approaches can perform such dimensionality reduction. For example, based on data 100 describing the operation of an HVAC system and the Boussinesq equation, the system performs a computational fluid dynamics (CFD) simulation 102 to generate an estimate of the thermal state at every point (or location) of the environment. In other words, the CFD simulation 102 technique can be used to simulate the dynamics of airflow in the environment. The thermal state at every point generated by the CFD simulation 102 is referred to as the full thermal state 104. The full thermal state 104 is considered as the ground truth value of the thermal state at every point of the environment. Furthermore, the ROM 106 is derived such that the estimate of the dynamics of airflow in the environment using the ROM 106 approximates the full thermal state 104.

[0036] For example, the CFD simulation 102 produces a complete thermal state 104 at each instant in time. The complete thermal state 104 at an instant in time is called a snapshot. The different snapshots of the CFD simulation 102 can be organized into a matrix that captures the airflow dynamics, and the ROM 106 is determined by solving an optimization problem such that the difference between the airflow dynamics estimated by the ROM 106 in a low-dimensional space and the airflow dynamics provided by the CFD simulation 102 is minimized. Such a scheme is useful for finding the minimum ||AX i -X i+1 It can be thought of as a linear estimation of ROM106 as a matrix A by solving the minimization problem given by ||, where A is ROM106, X i represents the thermal state at time i. The matrix A is used to linearly estimate the airflow dynamics in a low-dimensional space (e.g., x i+1 =A * x i), can be used. This minimization problem can be solved by several matrix decomposition techniques. For example, matrix A can be determined using proper orthogonal decomposition (POD) with Galerkin Projection, dynamic mode decomposition (DMD), spectral POD, and balanced POD, etc. Using the techniques described above, the minimization problem can be solved to determine matrix A, i.e., ROM106.

[0037] In some cases, the determined ROM 106 may not be linear. The determined ROM 106 may be a quadratic equation or a nonlinear ROM. Furthermore, the aforementioned techniques that can solve the minimization problem, alone or in combination, do not guarantee that the determined ROM 106 can be used for online HVAC control, for example, the determined ROM 106 may lack stability or observability. Because the ROM 106 operates in a low-dimensional space, an observer 108 is needed that can generate an estimate / reconstruction 112 of the thermal state at all desired points in the environment using the airflow dynamics captured by the ROM 106 in the low-dimensional space and actual measurements at several points in the environment. According to an embodiment, the observer 108 is a function of the ROM 106, i.e., the ODE and / or matrix A, and the placement of sensors in the environment. Specifically, the structure of the observer 108 includes a gain, referred to herein as the gain matrix K, and designing the observer 108 means determining the matrix K.

[0038] Some embodiments are based on the recognition that since matrix K depends on matrix A, matrix A is determined first, and then matrix K is determined based on matrix A. However, some embodiments are based on the recognition that the calculation of matrix K may affect the determination of matrix A. In other words, there is an interdependency between the determination of matrix A and the determination of matrix K. As a result, some embodiments are based on the recognition that a joint and interdependent optimization of matrix A and matrix K may improve the accuracy and global optimality of the thermal state reconstruction 112. To that end, the system performs a multivariate optimization including a joint optimization of the structure of the observer 108 (i.e., matrix K) and the structure of the ROM 106 (i.e., matrix A). That is, the ROM 106 and the observer 108 are determined together as optimization parameters of a cost function that uses the complete thermal state 104 as the ground truth.

[0039] Thus, the system optimizes the cost function to reduce the difference between the complete thermal state 104 of the environment, which is simulated based on the received data 100 and a physical model of the airflow including the PDE, and the thermal state 112 of the environment, which is reconstructed by the observer 108 using the ROM 106 from the values ​​of the complete thermal state at a set of locations. The optimization of the cost function is a multivariate optimization of the structure of the observer 108 and the structure of the ODE or ROM 106. According to an embodiment, the set of locations of the sensors for measuring the thermal state is predetermined.

[0040] Additionally or alternatively, some embodiments are based on the recognition that matrix K is a function of the placement of sensors in the environment. The sensors may measure one or a combination of temperature, velocity, and humidity of air in the environment. The placement of sensors can be represented as a binary matrix C with elements corresponding to locations in the environment. The elements of matrix C have a value of 1 in elements corresponding to the location of the sensors and a value of 0 otherwise. In some embodiments, matrix C, i.e., the locations of sensors 110, are predetermined by the designer of the HVAC system. However, in some other embodiments, matrix C may be unknown and / or may be a design parameter. Some embodiments are based on the recognition that matrix K depends on matrices A and C, and thus matrix C may be incorporated as an optimization parameter in a multivariable optimization. In other words, the set of the structure of the ODE or ROM 106, the structure of the observer 108, and the locations of sensors 110 are jointly optimized as part of a multivariable optimization, referred to as "co-design." The co-design approach not only results in an accurate reconstruction 112 of the thermal state, but also in optimal locations of sensors 110.

[0041] Some embodiments are based on another realization that when sensor placement is also part of the multivariable optimization, it is computationally efficient to perform the multivariable optimization via pairwise combinations such that the observer 108 is determined considering the ODE and the set of sensor 110 locations, the ODE is determined considering the observer 108 and the set of sensor 110 locations, and the set of sensor 110 locations is determined considering the observer 108 and the ODE. In this manner, the convergence of the multivariable optimization can be improved.

[0042] Additionally or alternatively, in some embodiments, a set of actuator positions of the HVAC system is known and represented by matrix B. In such a case, the set of actuator positions is not considered in the multivariable optimization. However, when the actuator positions are unknown, matrix B is incorporated as an optimization parameter in the multivariable optimization to determine the set of actuator positions of the HVAC system. To that end, in an embodiment, the multivariable optimization may correspond to a joint optimization of the structure of the ODE or ROM 106, the structure of the observer 108, the set of positions of the sensors 110, and the set of positions of the actuators.

[0043] Further, in some embodiments, the controller of the HVAC system can be modified with the optimized structure of the observer 108 and the optimized structure of the ODE. In particular, the reconstructed thermal state 112 is used for online HVAC control (e.g., closed-loop control) to achieve a desired state (e.g., thermal comfort for occupants or minimum energy consumption of the HVAC system).

[0044] For example, the average temperature or speed of the occupancy of a room, i.e., the portion of the room occupied by an occupant, can be considered as the desired state. At certain conditions of the HVAC system, the resulting average temperature and speed of the occupancy may differ from the desired state. Using an optimized structure of ODEs, the temperature and speed of the occupancy can be reconstructed, so that the difference between the current and desired values ​​can be measured. Any feedback controller can then be used to modify the output of the HVAC system, e.g., the HVAC speed and temperature, which may act as a controller, to get closer to what is sensed at the occupancy. Once the error (difference between desired and reconstructed) is evaluated, the control (HVAC speed, angle or temperature) can be adjusted to reduce this error. Even simple control strategies, such as PID to limit seeking, can be used to systematically reduce such errors.

[0045] FIG. 2A is a schematic diagram illustrating the conversion of the structure of the observer 108 and the structure of the ODE from a multivariable optimization 200 to a multivariable constrained optimization 202 according to some embodiments. In the multivariable optimization 200, the optimization of both the structure of the observer 108 and the structure of the ODE are the primary optimization objectives. In other words, in the multivariable optimization 200, the observer 108 and the ODE are given equal / balanced importance 204. However, some embodiments are based on the recognition that the importance of the ODE and the observer 108 in the objective of online HVAC control is not equal. In one embodiment, since the objective of the observer 108 is to enable the use of the ODE, the ODE is the primary objective and the observer 108 is a secondary objective, i.e., the observer 108 is subordinate to the ODE. Therefore, simply considering the ODE and the observer 108 jointly may not be optimal.

[0046] To alleviate such problems, multivariable optimization 200 is transformed into multivariable constrained optimization 202. Multivariable constrained optimization 202 provides an interdependent but unbalanced multivariable optimization of the ODE and the observer, where the ODE is the primary optimization objective and the observer 108 is the secondary optimization objective. In other words, in multivariable constrained optimization 202, the ODE is given higher importance compared to the observer 108. Thus, the ODE and the observer 108 are given unequal / unbalanced importance 206. To achieve the interdependence, the system jointly optimizes the structure of the observer 108 and the structure of the ODE, for example, as part of a single unit multivariable optimization. To achieve the unbalanced importance 206, some embodiments use multivariable constrained optimization in preference to the optimization of the ODE, where the optimization of the structure of the observer 108 acts as a constraint to the optimization of the ODE (208). In this manner, optimization of the observer 108 structure becomes subordinate to optimization of the ODE, since the observer 108 structure does not have an independent optimization objective, but serves only to constrain the structure of the ODE.

[0047] Furthermore, some embodiments are based on the recognition that when the multivariable optimization 200 is transformed into a multivariable constrained optimization 202 with an unbalanced optimization objective, the multivariable constrained optimization 202 can be used to address some potential deficiencies in the structure of the ODE and observer 108. For example, some approaches that sequentially determine the ODE and observer 108, i.e., first determine the ODE and then determine the observer 108 by considering the ODE, may make the matrix A of the ODE unstable. As used herein, stability is the property that a solution that starts near an equilibrium point will stay near the equilibrium point forever. For example, the solution of the ODE will always remain close to the steady-state solution when the HVAC system is slightly perturbed. Thus, the performance of the HVAC system depends on stability. Therefore, even if the PDE is unstable, there is still a need to determine a stable ROM for control / observation.

[0048] 2B is a schematic diagram illustrating constrained optimization for determining a stable ROM, according to some embodiments. Some embodiments are based on the recognition that since multivariate optimization adds flexibility to the optimization of the structure of the ODE 210, it may be advantageous to perform a multivariate optimization subject to a stability constraint 214 imposed on the structure of the ODE 210. Such optimization may be referred to as constrained optimization. Thus, constrained optimization may refer to a joint optimization of the structure of the ODE 210, imposed with the stability constraint 214, and the structure of the observer 212. Such constrained optimization results in a stable ROM. According to an embodiment, the stability constraint 214 may be a Lyapunov inequality.

[0049] Additionally or alternatively, some embodiments are based on the recognition that constraints imposed in constrained optimization may include constraints on the structure of the observer 212 without significantly increasing the computational complexity of the multivariate optimization. Examples of such constraints include observability constraints. As referred to herein, observability is a measure of how well the internal state of a system, e.g., the temperature of the environment, can be inferred from knowledge of its external outputs, e.g., the measurements of a thermostat. Enforcing observability on the structure of the observer 108 may therefore be advantageous, as it ensures that the temperature in the environment can be estimated with a given set of measurements.

[0050] 2C is a schematic diagram illustrating constrained optimization with stability constraints 214 and observability constraints 216 according to some embodiments. Here, the constrained optimization may correspond to a joint optimization of the structure of the ODE 210 with the stability constraints 214 and the structure of the observer 212 with the observability constraints 216. Examples of the observability constraints 216 include the Lyapunov equation for the observability Gramian with the observability matrix as the maximum rank, etc. Such constrained optimization results in a stable observable ROM.

[0051] However, some embodiments are based on the recognition that the conversion of ODE optimization to multivariate optimization, and in particular to multivariate constrained optimization, creates additional problems that are not present in the conventional single-variable optimization of the structure of ODEs. For example, the techniques described above that can solve the optimization of the structure of ODEs are not possible or practical for the multivariate constrained optimization and / or constrained optimization described above with reference to Figures 2A, 2B, and 2C.

[0052] According to some embodiments, the multivariate constrained optimization can be solved using semidefinite programming (SDP).

[0053] FIG. 3 is a schematic diagram for solving a multivariate constrained optimization based on SDP 300 according to some embodiments. SDP is a subdiscipline of convex optimization that concerns the optimization of a linear objective function on the intersection of a cone of a semidefinite matrix with an affine space. The system optimizes a cost function using SDP that estimates the structure of an observer and an ODE that are mutually dependent. Some embodiments are based on the realization that the multivariate optimization of the structure of an observer and the structure of an ODE is a convex problem. As a result, the multivariate constrained optimization 300 is also a convex problem. Thus, in some embodiments, the multivariate constrained optimization 300 may correspond to a constrained convex multivariate optimization 302.

[0054] Further, some embodiments are based on the recognition that SDP is a sub-discipline of convex optimization, and although multivariate optimization of observer-ODE structures is a convex problem, the dimensionality of the ODEs still imposes a potentially infeasible computational burden on SDP in convex domains. However, some embodiments are based on the recognition that large convex but constrained problems are harder to solve than large non-convex but unconstrained problems, because SDP may not scale well with the large dimensions of matrix A (even though the dimensions of matrix A are much smaller than those of CFD). To that end, some embodiments are based on the recognition that it is beneficial to convert the constrained convex multivariate optimization 302 to the more easily solvable unconstrained non-convex multivariate optimization 304.

[0055] According to some embodiments, an indirect SDP method is used to convert the constrained convex multivariable optimization 302 to an unconstrained non-convex multivariable optimization 304. In particular, the indirect SDP method uses a Cholesky decomposition to convert the constrained convex multivariable optimization 302 to an unconstrained non-convex multivariable optimization 304. The indirect SDP method is further used to solve the unconstrained non-convex multivariable optimization to optimize the structure of the observer 108 and the structure of the ODE. As a result, an optimized structure of the observer 108 and an optimized structure of the ODE are generated (306).

[0056] System Overview 4 is a block diagram illustrating a system for controlling the operation of an HVAC system arranged to condition an environment, according to some embodiments. The system 400 may have a number of interfaces connecting the system 400 with other systems and devices. For example, a network interface controller (NIC) 414 is adapted to connect the system 400 to a network 416 via a bus 412. Via the network 416, wirelessly or wired, the system 400 may receive data 100 indicative of the operation of the HVAC system. Additionally or alternatively, the data 100 indicative of the operation of the HVAC system may be received via an input interface 402.

[0057] The system 400 comprises a processor 404 configured to execute stored instructions and a memory 406 storing instructions executable by the processor 404. The processor 404 may be a single-core processor, a multi-core processor, a computing cluster, or any number of other configurations. The memory 406 may include a random access memory (RAM), a read-only memory (ROM), a flash memory, or any other suitable memory system. The processor 404 is connected to one or more input / output devices via a bus 412. Furthermore, the system 400 comprises a storage device 408 adapted to store different modules that store executable instructions for the processor 404. The storage device 408 may be implemented using a hard drive, an optical drive, a thumb drive, an array of drives, or any combination thereof.

[0058] The storage device 408 is configured to store a CFD simulation module 410a for generating the full thermal state 104. In other words, the CFD simulation module 410a is used to simulate airflow dynamics in the environment. The storage device 408 is configured to store an indirect SDP module 410b for solving the constrained convex multivariable optimization 302.

[0059] In some embodiments, the processor 404 is configured to optimize a cost function to reduce the difference between a complete thermal state 104 of the conditioned environment, simulated from the received data 100 based on a physical model of the airflow including PDEs, and a thermal state 112 of the conditioned environment, reconstructed by the observer 108 using a ROM 106 of the airflow including ordinary differential equations (ODEs) from values ​​of the complete thermal state at a set of locations. The optimization of the cost function is a multivariate optimization of the structure of the observer 108 and the structure of the ODE. The set of locations for measuring the thermal state is predetermined and provided to the multivariate optimization. The processor 404 is further configured to modify the controller 420 of the HVAC system with the optimized structure of the observer 108 and the optimized structure of the ODE. In some embodiments, the multivariate optimization is a constrained optimization of the structure of the ODE, subject to constraints imposed on the structure of the observer 108. In some other embodiments, the constrained optimization is subject to stability constraints imposed on the structure of the ODE.

[0060] The system 400 may further include an output interface 420. In some implementations, the system 400 is configured to submit the optimized structure of the observer 108 and the optimized structure of the ODE to the controller 420 via the output interface 420.

[0061] CFD Solutions

[0062]

number

[0063] At block 502, the processor 404 performs a discretization of the governing equations. In one embodiment, the governing equations are the turbulent transient Boussinesq equations, which are the incompressible Navier-Stokes equations coupled with the heat transfer equation, and are given by the following using Einstein notation:

[0064]

number

[0065] In a computational fluid dynamics (CFD) simulation, the governing equations (e.g., Navier-Stokes or energy equations) are first discretized. For example, in the finite volume method, the domain (e.g., a room containing furniture, appliances and occupants) is divided into various small computational entities called finite volumes. The number of such finite volumes (cells) is a function of the geometry, the speed of the HVAC, the size of the room, the desired accuracy, the method involving various physics such as turbulence or radiation, etc. In each cell, the PDE can be written as a linear algebraic equation whose solution is a function of the neighboring cells. Thus, to find the thermal state (velocity, temperature, etc.), which is the solution of the PDE to be solved in the CFD method, it is necessary to solve the equations for all the volumes. Since all the equations after discretization are linear, a linear solver is used at this stage, as shown in block 504.

[0066]

number

[0067]

number

[0068] Mathematical realization of codesign 6 is a schematic diagram illustrating a co-design for determining a stable observable ROM with the sensor's positions, according to some embodiments. A data-driven constrained optimization framework 600 is formulated for model reduction and stability of the identified model. The data-driven constrained optimization framework 600 includes aspects of a weighted least squares problem 602, a stability guarantee equation 604, and an observability guarantee 606.

[0069]

number

[0070]

number

[0071]

number

[0072]

number

[0073]

number

[0074]

number

[0075]

number

[0076]

number

[0077]

number

[0078]

number

[0079]

number

[0080]

number

[0081]

number

[0082]

number

[0083]

number

[0084]

number

[0085]

number

[0086]

number

[0087]

number

[0088]

number

[0089] Thus, in the limited memory BFGS scheme, the processor 404 computes the gradient of the cost function in block 700. Additionally, the processor 404 computes the changes in the design variables (A, C, K, etc.) in block 702. In block 704, the processor 404 computes the Hessian approximation, and in block 706, the processor 404 updates the design variables.

[0090]

number

[0091]

number

[0092] We need to choose a method for applying gradient descent by transforming the problem into a non-convex one with an analytical expression for the gradient. In this case, we use a quasi-Newton method, which has fast convergence.

[0093]

number

[0094]

number

[0095]

number

[0096]

number

[0097]

number

[0098] 8 is a flow chart illustrating a method for controlling operation of an HVAC system arranged to condition an environment, according to some embodiments. At block 800, the method includes receiving data indicative of operation of the HVAC system.

[0099] In block 802, the method includes optimizing a cost function to reduce the difference between a full thermal state of the conditioned environment (e.g., full thermal state 104) simulated from the received data based on a physical model of the airflow including partial differential equations (PDEs) and a conditioned thermal state (e.g., reconstructed thermal state 112) reconstructed by an observer (e.g., observer 108) using a reduced order model (ROM) of the airflow (e.g., ROM 106) including ordinary differential equations (ODEs) from the values ​​of the full thermal state at a set of locations. The optimization of the cost function is a multivariate optimization of the structure of the observer and the structure of the ODE. According to an embodiment, a set of locations for measuring the thermal state is predetermined and provided to the multivariate optimization.

[0100] At block 804, the method includes modifying a controller of the HVAC system using the optimized structure of the observer and the optimized structure of the ODE.

[0101] 9 is a diagram illustrating control of a heating, ventilation, and air conditioning (HVAC) system 910 using system 400, according to some embodiments. HVAC system 910 is positioned to condition a room 900. Room 902 is occupied by occupants 902, 904, 906, and 908. Arrow 914 represents air supplied by HVAC system 910 to condition room 900. System 200 may estimate thermal conditions at different locations in room 900. Specifically, in some embodiments, system 200 may estimate thermal conditions at the locations of occupants 902, 904, 906, and 908 within room 900.

[0102] In an embodiment, a controller 912 associated with the HVAC system 910 may be provided with target thermal conditions at the locations of the occupants 902, 904, 906, and 908. For example, a human operator may input the target thermal conditions at the locations of the occupants 902, 904, 906, and 908. Based on the thermal conditions at the locations of the occupants 902, 904, 906, and 908 estimated by the system 400, the controller 912 may generate control commands to achieve the target thermal conditions at the locations of the occupants 902, 904, 906, and 908. For example, the control commands may include a compressor speed of the HVAC system 910, a speed of an indoor fan of the HVAC system 910, and a position of an expansion valve of the HVAC system 910, etc. The controller 912 further controls the HVAC system 910 according to the control commands to achieve the target thermal conditions at the locations of the occupants 902, 904, 906, and 908.

[0103] The above description provides only exemplary embodiments and is not intended to limit the scope, applicability, or configuration of the present disclosure. Instead, the above description of exemplary embodiments provides those skilled in the art with an effective description for implementing one or more exemplary embodiments. Various changes may be made in the function and arrangement of elements without departing from the spirit and scope of the disclosed subject matter as set forth in the appended claims.

[0104] In the above description, specific details are provided to thoroughly understand the embodiments. However, it will be understood by those skilled in the art that the embodiments may be practiced without these specific details. For example, systems, processes, and other elements in the disclosed subject matter may be shown as components in block diagram form so as not to obscure the embodiments in unnecessary detail. In other instances, well-known processes, structures, and techniques may be shown without unnecessary detail so as not to obscure the embodiments. Additionally, like reference numbers and designations in the various drawings indicate like elements.

[0105] Also, each embodiment may be described as a process that is depicted as a flowchart, a flow diagram, a data flow diagram, a structure diagram, or a block diagram. Although a flowchart may describe operations as sequential processes, many of the operations may be performed in parallel or simultaneously. Also, the order of operations may be rearranged. A process terminates when its operations are completed, but may have additional steps not discussed or included in the diagram. Moreover, not all operations in a specifically described process may occur in all embodiments. A process may correspond to a method, a function, a procedure, a subroutine, a subprogram, or the like. When a process corresponds to a function, the end of the function may correspond to a return of the function to the calling function or to the main function.

[0106] Furthermore, embodiments of the disclosed subject matter may be implemented, at least in part, either manually or automatically. The manual or automatic implementation may be performed, or at least assisted, by the use of a machine, hardware, software, firmware, middleware, microcode, hardware description languages, or any combination thereof. When implemented in software, firmware, middleware, or microcode, the program code or code segments to perform the necessary tasks may be stored on a machine-readable medium. A processor(s) may perform the necessary tasks.

[0107] The various methods or processes outlined herein may be coded as software executable on one or more processors employing any one of a variety of operating systems or platforms. Furthermore, such software may be written using any of a number of suitable programming languages ​​and / or programming or scripting tools, and compiled as executable machine code or intermediate code that runs on a framework or virtual machine. Typically, the functionality of the program modules may be combined or distributed as desired in various embodiments.

[0108] The embodiments of the present disclosure may be embodied as a method, of which an example is provided. The acts performed as part of the method may be ordered in any suitable manner. Thus, while an exemplary embodiment shows acts as sequential, embodiments may be constructed in which acts are performed in a different order than illustrated, which may include performing some acts simultaneously. Although the present disclosure has been described with reference to certain preferred embodiments, it will be understood that various other adaptations and modifications can be made within the spirit and scope of the present disclosure. It is therefore the object of the appended claims to cover all such variations and modifications as fall within the true spirit and scope of the present disclosure.

Claims

1. 1. A system for controlling operation of a heating, ventilation and air conditioning (HVAC) system arranged to regulate an environment, comprising: at least one processor; and a memory, the memory, when executed by the at least one processor, providing the system with: receiving data indicative of the operation of the HVAC system controlled by a controller of the HVAC system; and optimizing a cost function to reduce a difference between a complete thermal state of the conditioned environment, which is simulated from the received data based on a physical model of airflow including partial differential equations (PDEs), and a thermal state of the conditioned environment, which is reconstructed by an observer from values ​​of the complete thermal state at a set of locations using a reduced order model of airflow (ROM) including ordinary differential equations (ODEs), wherein the optimization of the cost function is based on the set of locations for measuring the thermal state, the structure of the observer, and the structure of the ODEs. and a multivariate optimization for jointly and interdependently optimizing a structure of the ODE, the multivariate optimization being a constrained optimization of the structure of the ODE subject to constraints imposed on the structure of the observer, the constrained optimization being subject to stability constraints imposed on the structure of the ODE and subject to observability constraints imposed on the structure of the observer, the at least one processor optimizing the cost function using a semidefinite programming (SDP) that estimates the structure of the observer and the structure of the ODE that are interdependent, the system further comprising: storing instructions for modifying the controller of the HVAC system with the optimized structure of the observer and the optimized structure of the ODE, such that the operation of the HVAC system is controlled by the modified controller of the HVAC system based on the set of locations for measuring the thermal conditions; The multivariate optimization is a constrained convex problem, and the SDP uses an indirect solution to convert the constrained convex multivariate optimization to an unconstrained non-convex multivariate optimization, and solves the unconstrained non-convex multivariate optimization to optimize the structure of the observer and the structure of the ODE.

2. The system of claim 1 , wherein the SDP converts the multivariate optimization of the constrained convex problem to the unconstrained non-convex multivariate optimization using a Cholesky decomposition.

3. The system of claim 1 , wherein the cost function combines estimates of the observer, the ODE, and the set of positions.

4. 1. A method for controlling operation of a heating, ventilation and air conditioning (HVAC) system arranged to regulate an environment, comprising: receiving data indicative of the operation of the HVAC system controlled by a controller of the HVAC system; optimizing a cost function to reduce a difference between a complete thermal state of the conditioned environment, which is simulated from the received data based on a physical model of airflow including partial differential equations (PDEs), and a thermal state of the conditioned environment, which is reconstructed by an observer from values ​​of the complete thermal state at a set of locations using a reduced order model of airflow (ROM) including ordinary differential equations (ODEs); wherein the optimization of the cost function is a multivariate optimization that jointly and interdependently optimizes the set of locations for measuring the thermal state, a structure of the observer, and a structure of the ODE, the multivariate optimization being a constrained optimization of the structure of the ODE subject to constraints imposed on the structure of the observer, the constrained optimization being subject to stability constraints imposed on the structure of the ODE and subject to observability constraints imposed on the structure of the observer, and the cost function is optimized using a semidefinite programming (SDP) that estimates the structure of the observer and the structure of the ODE, which are interdependent; and the method further comprises: modifying a controller of the HVAC system with the optimized structure of the observer and the optimized structure of the ODE, such that the operation of the HVAC system is controlled by the modified controller of the HVAC system based on the set of locations for measuring the thermal conditions; The multivariate optimization is a constrained convex problem, and the SDP uses an indirect solution that converts the constrained convex multivariate optimization to an unconstrained non-convex multivariate optimization, and solves the unconstrained non-convex multivariate optimization to optimize the structure of the observer and the structure of the ODE.

5. The method of claim 4 , wherein the SDP transforms the constrained convex multivariate optimization into the unconstrained non-convex multivariate optimization using reduced dimensional approximations.

6. The method of claim 5 , wherein the reduced dimensional approximation comprises a Cholesky decomposition.

Images