Improvements in control systems for thermal conditions

Abstract

A method for offline Koopman system identification is provided which includes generating snapshots representing a thermal comfort quantification index, reshaping these snapshots into 3D tensors, and splitting them into training, validation, and testing datasets. An encoder is used to reduce state space shape, while a dynamic auxiliary network calculates the A and B matrices representing a reduced order system. The auxiliary output advances the reduced order system into a time step with the control of the boundary condition. A decoder is used to reconstruct the system to its original state for another time step, and a loss function for reconstruction, linearisation and prediction is computed.

Description

IMPROVEMENTS IN CONTROL SYSTEMS FOR THERMAL CONDITIONS Technical Field

[0001] The present invention relates to improvements in feedback control systems for thermal conditions and in particular heating, ventilation, and air conditioning (HVAC) feedback control systems. Background of Invention

[0002] Control and monitoring for indoor thermal comfort is desirable within indoor spaces such as office buildings as well striking a balance between comfort and energy consumption. Traditional feedback control systems for thermal conditions rely on temperature readings from thermometers as an indication of the indoor thermal comfort – however it is desirable to define thermal comfort from the point of view of an occupant within a space. That is, thermal comfort is achieved when an occupant’s body within the space reaches equilibrium with the surrounding conditions. In that regard, indices have been developed which define thermal comfort conditions, notably Predicted Mean Vote (PMV) which attempts to account for the physiological conditions of occupants such as clothing insulation and metabolic rate and air characteristics such as air radiant and dry bulb temperatures, relative humidity, air velocity, and air pressure.

[0003] Systems which utilise a conventional PMV index model tend to be highly complex (since the model includes nonlinear recursive relations and requires iterative calculation of the roots of nonlinear equations for example).

[0004] “Black Box” optimizations have also been proposed for approximating and predicting thermal comfort conditions by integrating backpropagation of neural networks to ensure that the system is adapted to user’s thermal preferences. However, these arrangements have difficulty in achieving real-time results and have delays in tracking for indoor conditions.

[0005] Attempts have been made at a convergence between a PMV index, represented by using Computational Fluid Dynamics and embedded simple feedback control models but the coupling of control and constant system modelling requires extensive computational resources.

[0006] Another problem is that feedback control systems for thermal conditions are inherently nonlinear - which is computationally intensive. Attempts have been made to “linearise” feedback control systems for thermal conditions using model linearisation techniques, for example Generalised Laplace Analysis (GLA), Proper Orthogonal Decomposition (POD)-Galerkin and Dynamic Mode Decomposition (DMD), but these model reduction methods tend to be applied to specific airflow characterisation cases with limited application to ventilation systems.

[0007] It would therefore be desirable to provide a linear model of the nonlinear thermal comfort index. It would further be desirable to provide a method and system which ameliorates or at least alleviates one or more of the above-mentioned problems or provides a useful alternative.

[0008] A reference herein to a patent document or other matter which is given as prior art is not to be taken as an admission that that document or matter was known or that the information it contains was part of the common general knowledge as at the priority date of any of the claims. Summary of Invention

[0009] According to a first aspect, the present invention provides, a method for offline Koopman system identification, including the steps of: a) providing one or more generated snapshots ( ^^) representative of a thermal comfort quantification index, each being associated with a fixed range of boundary conditions for at least one of temperature, velocity ( ^^) or pressure (p); b) reshaping the one or more generated snapshots ( ^^) to generate one or more 3D tensors having variables (No. trajectory, ^^^^, ^^ ) and (No. trajectory, ^^^^, ^^) where: No. trajectory is the trajectory number; ^^^^is Timesteps for CFD horizon, sensor data time range or thermal images’ timesteps; ( ^^) are generated snapshots; and ( ^^) is a range of boundary conditions for temperature and velocity; c) splitting the generated 3D tensors into at least one training dataset, at least one validation dataset and at least one testing dataset; d) providing an encoder to reduce state space shape from ^^ to ^^^^= ^^^^(^^^^); e) providing a dynamic auxiliary network having one or more layers to thereby calculate an A matrix and a B matrix, the A and B matrix being representative of a reduced order system definition; f) computing the auxiliary output of ^^̂^^+ ^^= ^^ ^^^^+ ^^ ^^ to … advance the reduced order system into atime step with the control of the boundary condition; g) providing a decoder such that the reduced order system may be reconstructed to the original state for a further time step from ^^̂^^+ ^^→ ^^^^+ ^^; and h) computing a loss function for reconstruction, linearisation and prediction.

[0010] Preferably, the method further includes training steps, wherein, over a predetermined number of epochs, for each of steps c) to h), the method further includes: i) calculating a loss function for the training dataset; j) calculating a loss function for the validation dataset; k) optimising the training process via an optimizer; and l) passing through stopping conditions of training loss vs validation loss.

[0011] Preferably, the method further includes the steps of: m) saving encoder hyperparameters A and B matrices; n) testing prediction accuracy against testing dataset for the autoencoder.

[0012] Preferably, the generated snapshots ( ^^) that are representative of a thermal comfort quantification index include Predictive Mean Vote - Computational Fluid Dynamic (PMV-CFD).

[0013] The range of boundary conditions may further include relative humidity.

[0014] Preferably, at step b), reshaping includes transforming 3D meshes to a point cloud representation, thereby reducing computational resources.

[0015] Preferably, at step e), the dynamic auxiliary layer is created by way of:and ^^^^and ^^^^^^are orthogonal unitary columns and ∑ ^^ ^^= ^^ ^^ ^^ ^^( ^^^^, ^^)^^=1is a diagonal hierarchical matrix.

[0016] In an alternative, at step e), the dynamic auxiliary layer is created by way of:where: ^^^^are eigenvalues represented by complex imaginary and real values ^^±= ^^ ± ^^ ^^ , denoting the oscillating and decaying nature of airflow and are parameterised by the cosine and sine matrix along with exponential representation.

[0017] Preferably, at step k), the optimizer includes one or more of an ADAM optimizer, Adadelta optimiser or Stochastic Gradient Descent (SGD) optimiser.

[0018] According to a second aspect, the present invention provides a method for generating an online feedback controller, including the steps of: a) embedding an encoder for ^^^^( ^^^^ ^^ ^^) → ^^^^ ^^ ^^, where: ^^^^, is a trained encoder network with optimal hyperparameters; b) applying a Linear-Quadratic Regulator (LQR) function on a predetermined set of Koopman operator matrices and a quadratic and regulator cost matrices to provide a K gain matrix; and c) setting the system environment to an initial condition ^^^^.

[0019] Preferably, at step a), the trained encoder network with optimal hyperparameters is the output of the offline Koopman system identification from the first aspect of the present invention.

[0020] Preferably, at step b), the predetermined set of Koopman operator matrices are the matrices A and B, derived from the first aspect of the present invention.

[0021] Preferably, at step b), the K gain matrix is provided by the function ^^ ^^ ^^( ^^, ^^, ^^, ^^) → ^^ ^^ ^^ ^^. In an alternative, at step b), the K gain matrix is provided by an Model Predictive Control (MPC) function.

[0022] Preferably, for each t in timesteps for control ^^^^, the method further includes the steps of: d) embedding an encoder for the current state ^^^^( ^^) → ^^; e) applying optimal feedback control ^^ ^^ ^^ ^^( ^^ − ^^^^ ^^ ^^); and f) adding disturbances updates online. Brief Description of Drawings

[0023] The invention will now be described in further detail by reference to the accompanying drawings. It is to be understood that the particularity of the drawings does not supersede the generality of the preceding description of the invention.

[0024] Figure 1 is a flow diagram illustrating the process steps for a method for offline Koopman system identification in accordance with an exemplary embodiment of the present invention;

[0025] Figure 2 is a flow diagram illustrating the process steps for a method of generating an online feedback controller in accordance with an exemplary embodiment of the present invention;

[0026] Figure 3 is a system diagram illustrating the operation of the methods described in Figure 1 and Figure 2;

[0027] Figure 4 is a schematic diagram illustrating universal system identification of Koopman operators with a variety of observable behaviour constituting the dictionary of the trajectories in accordance with an exemplary embodiment of the present invention;

[0028] Figure 5 is a schematic diagram illustrating an example of a Koopman operator integrated with an autoencoder deep network in accordance with an exemplary embodiment of the present invention;

[0029] Figure 6 is a schematic diagram illustrating an example of 3D Geometry of an office space;

[0030] Figure 7 is a schematic diagram illustrating PMV index variables and ISO7730 categories in accordance with an exemplary embodiment of the present invention;

[0031] Figures 8a, 8b, 8c and 8d are point cloud diagrams illustrating PMV index point clouds for transient (a) ^^^^= 1, (b) ^^^^= 50, (c) ^^^^= 200 and (d) ^^^^= 600;

[0032] Figure 9 illustrates charts for autoencoder hyperparameter tuning in accordance with an exemplary embodiment of the present invention;

[0033] Figure 10 is a chart illustrating network training versus validation performance through defined epochs in accordance with an exemplary embodiment of the present invention;

[0034] Figure 11 is a chart illustrating eigenvalue visualisation for the real plane in accordance with an exemplary embodiment of the present invention;

[0035] Figure 12a is a PMV-CFD cloud point diagram illustrating the second time step in a testing trajectory dataset in accordance with an exemplary embodiment of the present invention;

[0036] Figure 12b is a cloud point diagram showing Koopman autoencoder network cloud point predictions for an exact time step in accordance with an exemplary embodiment of the present invention;

[0037] Figure 13 is a chart illustrating PMV index prediction error relative to the CFD global state results with the Koopman parameterisation method and SVD method in accordance with an exemplary embodiment of the present invention;

[0038] Figure 14 is a chart showing network with control input training vs validation losses for 200 epochs in accordance with an exemplary embodiment of the present invention;

[0039] Figure 15 is a schematic diagram of a closed loop latent linear system control;

[0040] Figure 16 is a chart illustrating closed loop performance on the latent space points for the control horizon;

[0041] Figure 17A to 17D are charts illustrating CFD-PMV points tracking with the application of various Koopman control inputs in accordance with an exemplary embodiment of the present invention;

[0042] Figure 18A and 18B is a schematic diagram of a sample of CFD tracked points’ locations given in Figure 17A to 17D and illustrates Transitioning between 2D CFD-PMV plane distribution and sample points of CFD point cloud distribution in accordance with an exemplary embodiment of the present invention; and

[0043] Figure 19 is a chart illustrating closed-loop performance of PMV-RC lumped LQR controller in accordance with an exemplary embodiment of the present invention. Detailed Description

[0044] The present invention may be utilised in the context of thermal comfort within a building, and it will be convenient to describe the invention in relation to that exemplary, but non-limiting, application. It will be appreciated that the present invention is not limited to that application and may for example, be applied beyond buildings and indeed to any confined environment including electric vehicles, aviation and other critically controlled environments. Other applications include indoor air quality characterisation with nonlinear distribution of particles and gas concentration. In addition, the invention may utilise several data driven approaches as basis of the neural network input, such as thermal images provided by thermal cameras and optimally located sensors, for example.

[0045] The present invention allows for a linear model of the thermal comfort index (which is strongly nonlinear) by way of Koopman eigenfunctions (which remain interpretable in high dimensional systems) and deep learning which avoids the need for a full representation of the system by way of CFD. In the present invention, thermal comfort dynamics may be modelled by: (1) the intrinsic coordinates of the thermal comfort distribution being defined by an autoencoder neural network for reconstruction of the dynamics, and (2) detecting and predicting the future thermal comfort state. The validity of the autoencoder prediction may be compared to exact CFD simulation results with a relative error calculated. Advantageously, the present invention allows for a more accurate offline estimation of the thermal comfort model compared to the constant update of an online black box method or by way of pure CFD modelling (which is computationally expensive).

[0046] Referring to Figure 1, there is shown a flow diagram illustrating the process steps for a method 100 for offline Koopman system identification, while Figure 2 shows a flow diagram illustrating the process steps for method 200 of generating an online feedback controller. The methods 100, 200 may be run via one or more servers [not shown] having one or more databases [not shown] connected to one or more computing devices [not shown] and may be communicatively coupled to a cloud computing environment and interconnected via a network such as the internet or a mobile communications network.

[0047] Although “cloud” has many connotations, according to embodiments described herein, the term includes a set of network services that are capable of beingused remotely over a network, and the method described herein may be implemented as a set of instructions stored in a memory and executed by a cloud computing platform. The software application may provide a service to one or more servers or support other software applications provided by a third-party servers. Examples of services include a website, a database, software as a service, or other web services. Computing devices may include smartphones, tablets, laptop computers, desktop computers, server computers, among other forms of computer systems.

[0048] The transfer of information and / or data over the network can be achieved using wired communications means or wireless communications means. Also, embodiments need not take place over a network, and the method steps could occur entirely on a client or server processing system.

[0049] At step 105, the method 100 begins where one or more generated snapshots ( ^^) representative of a thermal comfort quantification index, each being associated with a fixed range of boundary conditions for at least one of temperature and velocity ( ^^) are provided.

[0050] Control then moves to step 110, where the one or more generated snapshots ( ^^) are reshaped to generate one or more 3D tensors having variables (No. trajectory, ^^^^, ^^ ) and (No. trajectory, ^^^^, ^^), where: No. trajectory is the trajectory number, ^^^^is Timesteps for CFD horizon, ( ^^) are generated snapshots and ( ^^) is a range of boundary conditions for temperature and velocity. The “No. trajectory” is a hyperparameter that is suggested and later optimised during hyperparameter tuning. This is calculated as CFD time horizon= No.trajectory * Ts (tensor’s time steps).

[0051] Control then moves to step 115, where the generated 3D tensors are split into at least one training dataset, at least one validation dataset and at least one testing dataset. Splitting of the generated 3D tensors may be carried out in any suitable manner, for example, a Scikit Learn train and test model split (to split the dataset arrays into random train and test datasets based on percentages). These percentages may include 60% for training, 20% for validation and 20% testing for example.

[0052] An encoder is then provided at step 120, to reduce state space shape from ^^ to ^^^^= ^^^^(^^^^);

[0053] At step 125, a dynamic auxiliary network is provided, having one or more layers, to thereby calculate an A matrix and a B matrix, the A and B matrix being representative of a reduced order system definition.

[0054] Control then moves to step 130 where the auxiliary output of ^^̂^^+ ^^= ^^ ^^^^+ ^^ ^^ is computed to advance the reduced order system into a time step with the control of the boundary condition.

[0055] A decoder is provided at step 135 such that the reduced order system may be reconstructed to the original state for a further time step from ^^̂^^+ ^^→ ^^̂^^+ ^^. At step 140, a loss function for reconstruction, linearisation and prediction is computed.

[0056] The method preferably further includes training steps, wherein, over a predetermined number of epochs, for each of steps 115 to 140, the method further includes: calculating a loss function for the training dataset; calculating a loss function for the validation dataset; optimising the training process via an optimizer; and passing through stopping conditions of training loss versus validation loss. The optimizer may include one or more of an ADAM optimizer, Adadelta optimiser or SGD optimiser.

[0057] Preferably, the method the method further includes the steps of: saving encoder hyperparameters A and B matrices and testing prediction accuracy against testing dataset for the autoencoder.

[0058] It will be appreciated at step 105 the generated snapshots ( ^^) that are representative of a thermal comfort quantification index include Predictive Mean Vote - Computational Fluid Dynamic (PMV-CFD). Further, at step 105, the range of boundary conditions may further include relative humidity or generated thermal images.

[0059] At step 110, reshaping includes transforming 3D meshes to a point cloud representation, thereby reducing computational resources.

[0060] At step 125, the dynamic auxiliary layer may be created by way of:where:and ^^ and ^^^^are orthogonal unitary columns and ∑ = ^^ ^ ^^ ^^ ^^ ^^^ ^^ ^^( ^^^^, ^^)^^=1is a diagonal hierarchical matrix. In an alternative, at step 125, the dynamic auxiliary layer may be created by way of:where the eigenvalues ^^^^would be represented by complex imaginary and real values ^^±= ^^ ± ^^ ^^ , denoting the oscillating and decaying nature of the airflow streamlines for the indoor zone. Typically, these oscillations and damping are parameterised by cosine and sine matrix along with exponential representation, allowing the eigenvalues to vary across the transient state.

[0061] Figure 2 is a flow diagram of the process steps for method 200 for generating an online feedback controller in accordance with an exemplary embodiment of the present invention. At step 205, an encoder is embedded for ^^^^( ^^^^ ^^ ^^) → ^^^^ ^^ ^^, where: ^^^^, is a trained encoder network with optimal hyperparameters.

[0062] Control then moves to step 210 where a Linear-Quadratic Regulator (LQR) function is applied on a predetermined set of Koopman operator matrices and quadratic and regulator cost matrices to provide a K gain matrix. Control then moves to step 215, where the system environment is set to an initial condition ^^^^. Preferably, backward propagation is conducted during the training process to optimise the weights and biases definition.

[0063] Preferably, the encoder is applied on a suggested matrix for the desired state of the PMV index for each point for instance PMV index is preferred to reach 0 as much as possible to achieve optimal thermal comfort. The suggested ^^^^ ^^ ^^is the desired PMV values for each point with the same size as the original ^^. In this case, driving the original state of PMV as close as possible to the desired state of PMV.

[0064] Preferably, at step 205, the offline system identification calculates the A and B matrices required for the LQR application.

[0065] Preferably, at step 205, the trained encoder network with optimal hyperparameters is the output of the offline Koopman system identification obtained from the method of Figure 1. The hyperparameters may be determined initially after running the network for the first time. The tuning of hyperparameters may be carried out via any suitable arrangement (e.g. Ray Tune or the like), which enables suggesting a range or a matrix of the hyperparameters of the network and runs the network to test which set gives the best training loss faster. The hyperparameters may include batch size, number of trajectory, learning rate value and loss hyperparameters.

[0066] After training the network and optimising the hyperparameters, the saved network is tested with the testing dataset similar to the training process with the loss calculations to obtain a prediction and to calculate the relative error between this exact prediction time step and the CFD generated timestep.

[0067] Preferably, at step 210, the predetermined set of Koopman operator matrices are the matrices A and B, derived from the method of Figure 1.

[0068] Preferably, at step 210, the K gain matrix is provided by the function ^^ ^^ ^^( ^^, ^^, ^^, ^^) → ^^ ^^ ^^ ^^. After the definition of the system A and B matrices, a Q cost matrix penalises the state space to reach the desired state space while the R matrix penalises the energy consumption of the control input.

[0069] For each t in timesteps for control ^^^^, the method of Figure 2 further includes the steps of: embedding an encoder for the current state→ ^^; applying optimal feedback control ^^ ^^ ^^ ^^( ^^ − ^^^^ ^^ ^^); and adding disturbances updates online.

[0070] It will be appreciated that LQR may be employed for testing the suitability of the Koopman operator linearisation approach for optimising thermal conditions for an indoor environment – for example its capacity to achieve a balance between desired PMV conditions and minimising control effort. In addition, LQR optimal control problem tracks the state space to the reference / desired state. The LQR control is primarily designed to regulate the behaviour of the zone inlet of the supply system with a quadratic loss function. The quadratic loss function may be given as:where, ^^ is a positive semidefinite matrix that weighs the deviation of the state space, and ^^ is a positive definite matrix that weighs the control effort. Based on that equation, the system may be capable of recovering previous time steps and to predict the next state spaces. In other words, the system should be tractable, thus, applying the optimal control directly is taxing. With the application of the Koopman operator autoencoder, the nonlinearity of the PMV input may be mapped to linear dynamics in the latent space. Therefore, the application of optimal control in the auxiliary control network is ^^( ^^, ^^). Conformed to the Koopman linearisation model and given the desired state with the ISO 773 PMV (for example) satisfactory categories as ^^^^^^+^^1^^, the LQR optimal control problem may be rewritten as:Hence, the optimal LQR control input may take the form ofIn the context of solving the optimal control LQR problem, Pontryagin’s Minimum principle may be used, which states that the optimal control minimises a Hamiltonian function and leads to optimality. To solve the Hamiltonian function, discrete-time Ricatti equation provides the solution for the infinite horizon time-invariant aswhere ^^ is a symmetric matrix, called the state feedback gain matrix with ^^ and ^^ as described above and once the ^^ is found from the Ricatti equation, the control input in can be rewritten by induction backward as

[0071] This control policy may then be applied to the PMV-CFD based model previously described.

[0072] Figure 3 is a system diagram illustrating the operation of the methods described in Figure 1 and Figure 2. Shown is a system 300 having an offline Koopman system identification component 305 and an online feedback controller component 325. Sensor data buffer storage 305 includes either historical data or thermal PMV index data (which may be collected from optimally located sensors, stored CFD snapshots, or thermal images). This information provides input for the Koopman identificationprocess 310 where it is reshaped as previously mentioned with reference to Figure 1. The Koopman identification process 310 processes the input data by applying the linearisation and model reduction of Koopman operator and SVD auxiliary layer and calculates A and B matrices. The A and B matrices are input into the controller design 320 with the estimation of Q and R cost matrices. This process may be conducted offline to calculate the K gain matrix required for LQR feedback controller or MPC optimisation input.310, trained offline, contributes with tuned hyperparameters to the online encoder stage 330. Moving to online feedback controller component, system snapshots 345 provide input ^^ to the online encoder 330, utilising tuned hyperparameters from 310. The encoded states ^^ are incorporated with the output of 320 in order to reach the desired states considering the cost matrices and A and B matrices through ^^ ^^ ^^ ^^( ^^ − ^^^^ ^^ ^^) of 335. The calculated control signal is combined with consideration of disturbances 340 which may be presented as a change in occupancy number, metabolic rate change and / or humidity change, for example. The optimal control signal is applied to the system to improve the thermal conditions.

[0073] Figure 4 is a schematic diagram 400 illustrating universal system identification of Koopman operators with a variety of observable behaviour constituting the dictionary of the trajectories in accordance with an exemplary embodiment of the present invention. The Koopman operator parameterisation shown in Figure 4 is used as the auxiliary nonlinear mapping and added to the autoencoder network as will be described with reference to Figures 9 to 13.

[0074] The Koopman operator allows for full nonlinear representation of high dimensional dynamics being captured by a linear setting but without sacrificing required information which would otherwise compromise the performance of the model.

[0075] The Koopman operator, which is the span of the eigenfunctions associated with eigenvalues near a unit circle in discrete time or near the imaginary axis in continuous time, is utilised in the present invention. These eigenfunctions describe the long-term dynamics of observables that emerge after fast transients. When a spectral gap exists, this has subsided and can be used in low-dimensional approximations of the otherwise infinite-dimensional operator, which delineates fast and slow temporal dynamics. A set of Koopman modes is capable of reconstructing the state of systemsas a linear combination of the Koopman eigenfunctions with vectors - allowing for global linearisation to the original nonlinear state space.

[0076] In the present invention, Koopman analysis is utilised, which generalizes linear mode analysis from linear systems to nonlinear systems while preserving global nonlinear features.

[0077] An example of Koopman operator eigenfunctions and mode identification is provided below – where a dynamic system evolving on the state space of ^^ presented with the form of state input and output equations as: ^^^^+1= ^^( ^^^^), (1) where ^^ is the nonlinear dynamics that map the system state forward in time with ^^ as an integer index representing the instance of the state, and ^^ is the state of the system.

[0078] ^^ is not identified for the thermal comfort model development and the transient PMV measurements are available in representing ^^^^in Equation (1). This arrangement aims to determine spatiotemporal nonlinear pattern formation of the thermal comfort system with the Koopman operator. The Koopman operator is introduced as a linear operator ^^ that acts on scalar-valued functions of observables ^^ on the state space ^^. The concrete interpretation of the observables ^^ is the PMV CFD data which is mapped as the following ^^: ^^ → ℝ (as shown as 405a, 405b….405n). In theory, ^^ belongs to the state space where the Koopman operator describes the trajectory ( ^^1, ^^2, ...) and the span of the dictionary matrix contains all the observables as ^^ =(as shown as 410) ^^ maps ^^ into a space given by composition of the observable ^^ and the iterated map of ^^ as ^^ ^^(^^^^)= ^^(^^^^+1)(2)

[0079] ^^ (as shown as 415) is considered as a linear operator because it is homomorphic; addition and scaling are preserved since ^^(^^ ^^1+ ^^ ^^2)(^^)= ^^ ^^ ^^1(^^)+ ^^ ^^ ^^2( ^^) for any observable functions ^^1, ^^2∈ ^^, And scalars ^^, ^^; however, ^^ is an infinite dimensional operator, with a nonlinear dynamic system progressing on a finitedimensional space ^^. The Koopman operator ^^ can be extended to the large finite dimensional observables by using the eigenfunctions and eigenvalues of ^^, such that ^^^^: ^^ → ℝ represents the eigenvectors and ^^^^: ^^ → ℂ represents the eigenvalues computed through ^^ − ^^^^Ι = 0 where Ι is denoting the identity matrix (1) as ^^ ^^^^= ^^^^^^^^(3) ^^ = ^^^^^^^^^^^^−1(4)

[0080] Hence the vector-valued of observables ^^ in Equation (2) can be expressed in terms of a linear combination eigenvectors ^^^^and the eigenvalues ^^^^of ^^ as

[0081] In addition, ^^ can be expressed in terms of eigenfunctions ^^^^and eigenvalues: ^^^^( ^^^^+1) = ^^ ^^^^( ^^^^) = ^^^^^^^^( ^^^^) (7)

[0082] Leading to represent the observables ^^(^^^^+1)from Equation (5) as

[0083] This representation is advantageous as the dynamics associated with each eigenfunction are determined by its corresponding eigenvalue. Determining the eigenfunctions ensures to span an invariant subspace, subsequently restricting the Koopman infinite dimensionality to this subspace. By writing ^^ in terms of eigenfunctions, the complex evolution of(^^)can be substituted with linear evolution ^^^^, as a result, one can either solve a finite dimensional but nonlinear problem or an infinite dimensional but linear problem when eigenvalue, eigenfunctions and modes are computed. Preferably, the dictionary of observables is representative of the contributing dynamics of the thermal comfort quantification, hence, the performance of the present invention relies on assessing and evaluating the observables dictionary in 420 and testing the predictions of 425 of the estimated Koopman operator in line with the testing phase of the autoencoder performance.

[0084] For the quadratic nonlinear system of the PMV index, the eigenvalues ^^^^are be represented by complex imaginary and real values ^^±= ^^ ± ^^ ^^ , denoting the oscillating and decaying nature of the airflow streamlines for an indoor zone such as an office space or the like.

[0085] Typically, these oscillations and damping are parameterised by cosine and sine matrix along with exponential representation following Equation (8) allowing the eigenvalues to vary across the transient state.Singular Value Decomposition (SVD)

[0087] To compare the performance of the Koopman operator parameterisation and performance when dimensionally reduced, Principal Component Analysis (PCA) with truncated Singular Value Decomposition (SVD) is applied to further reduce the cost of computation and create a compressed version of the encoder output. This is applied to define the bases for the domain and range of the linear map associated with matrix multiplication. The training instance is the output of the encoder network ^^(^^) ∗∗and the associated reduced states ^^̂. These states are decomposed to into ^^^^̂components ^^̂^^, ^^, ^^ = 1, … , ^^^^̂. The SVD is computed on these states as the following:

[0088] ^^ and ^^^^are orthogona ∑ ^^ ^^ ^^l unitary columns and^^= ^^ ^^ ^^ ^^( ^^^^, ^^)^^=1is a diagonal hierarchical matrix representing the importance (i.e. the contribution of the state intrinsic coordinates to the nonlinearity of the entire system) of each state space. The columns of ^^^^represent the direction onto which the data is projected with the least approximation error. In addition, the first column of ^^^^^^1, ^^denotes the reducedrepresentations best expresses the variance in the training encoded data as Equation (12).

[0089] The amount of this variance is related to the ^^^^, ^^and the quantity can serve as a proxy for how much the variance is represented by the ^^ basis vectors of ^^^^.The smaller the number of basis vectors is the larger amount of training variance can be captured. This quantity can be determined by Equation (13).

[0090] In this arrangement, lower dimensional space can be determined in further detail as will be described further with reference to Figure 8.

[0091] Figure 5 is a schematic diagram 500 illustrating an example of a Koopman operator integrated with autoencoder deep network in accordance with an exemplary embodiment of the present invention.

[0092] The present invention utilises deep neural networks having the ability to efficiently represent and learn by itself complicated functions learned from data. The autoencoder is a self-supervised deep neural network that has a symmetric layer structure with a bottleneck in the middle, with an encoder as the first half and decoder the second half, as shown in Figure 5. At 505, the PMV snapshots are used as input values through training, validation and / or testing datasets through the encoder linear layers’ proposed structure. The reduced states ^^^^are then analysed and decomposed using both Koopman parameterisation and SVD auxiliary layer in 510, which enables the advancement of the reduced states through time to the next timesteps in ^^^^+1, as ^^^^+1= ^^ ^^^^+ ^^ ^^̂^^.The A matrix related to the state spaces are identified as the weights of the network proposed by the training process. In addition, the control input of the boundary conditions is introduced in order to calculate the B matrix. Next, ^^^^+1is reconstructed through the linear layers of the decoder to ^^^^+1the original size of the next timesteps as predictions as shown at 515. After imposing the loss functionsdescribed below at 520 the predictions are compared to the exact ^^^^+1from the CFD readings at 525 to ensure the accuracy of the predictions through iterative training. After reducing the training losses to a specified threshold, the system identification is complete at 535 through the identification of Koopman operator A and B matrices. Autoencoders generalize the idea of encoder and decoder beyond the deterministic functions to stochastic mapping. In the present invention the autoencoder network is leveraged to discover the Koopman operator eigenfunctions, while avoiding overfitting, addressing continuous spectrum and reducing the nonlinearity of the PMV index behaviour. With these three requirements, the network corresponds to three phases during training to reconstruct its input with minimal regression errors constrained with loss calculations and back propagation of the neural network features: 1. Intrinsic coordinates ^^^^definition for the input system ^^^^reconstruction spanned by the Koopman eigenfunctions ^^^^from Equation (7) as ^^^^= ^^^^(^^^^)as the encoder, and the inverse ^^ = ^^^^−1(^^^^)as the decoder to ensure system recovery. For this phase, the loss calculation as the mean square error of2. Koopman operator dynamics parametrisation with eigenvalues ^^^^from Equation (8) learned by an auxiliary network following Koopman parameterisation model or Singular Value Decomposition technique. The linearity loss is calculated as ^^^^=3. Future state prediction with the linear dynamics through the intrinsic coordinates, with the last term of neural network losses calculated as ^^^^ ^^ ^^ ^^=

[0093] The LQR controller is introduced through penalising the states ^^^^to reach the desired states ^^^^^^ ^^ ^^with the consideration of the cost matrices of the LQR controller as+ ^^̂^^^^^^̂ ^^̂^^. Hence the optimal control takes the form of ^^̂^∗^= ^^^^ ^^ ^^( ^^^^− ^^^^^^ ^^ ^^), which can be solved by a discrete Ricatti equation for example. Similarly, MPC controller can be introduced with the set points definition and optimisation.

[0094] To train the network, the data generated from the transient PMV CFD model (further described with reference to Figure 6) forms the training representing 60% of the PMV CFD data, validation and testing sets, each representing 20% of the PMV CFD transient data, the time instances for each dataset are not repeated following a distribution for sampling. The network is developed on the training set and compared to the validation set to avoid overfitting and tested for accuracy on the testing set as will be described with reference to Figure 8.

[0095] Figure 6 is a schematic diagram 600 illustrating an example of 3D Geometry of an office space. Thermal comfort and offline system identification is required to complete the loop of either efficient control of temperature or optimal sensor placement. Predicted Mean Vote index is established by a 3D CFD model, using, for example, incompressible Navier Stokes equations for the conservation of mass, momentum and energy. This model supplies the required observables to constitute the input dictionary of the deep neural network. To consider an accurate representation of the thermal comfort vertical and horizontal distribution, a simulation solver may be utilised (e.g. ANSYS R22.2 FLUENT for example) to develop the solution of the associated partial differential equations. In the example of Figure 6, the 3D model’s fixed geometry volume is a single occupied office 605 having dimensions of 3 meters in length × 3 meters in width × 2.6 meters in height, with one external wall 610 as a stationary wall of 30 degrees Celsius, to represent summer conditions, and interior walls 615 of 25 degrees Celsius, with the ceiling 620 and floor 625 considered as adiabatic surfaces. A Computer Simulated Person (CSP) 630 is preferably used to ensure the accuracy of the PMV index around a near-occupant region, and recognising the thermal plume of the occupant with heat generation of 63.96 W / m2. A Multizone polyhedral mesh is used to discretise the airflow domain to capture the irregular geometry of the occupant. A polyhedral mesh is used to reduce the number of iterations required to reach the CFD model convergence while using less elements when compared to a tetrahedral mesh and capturing high accuracy of the CSP geometry. A visualisation of the mesh is shown in Figure 6. The domain is divided into ^^^^ ^^= 1,048,576 non-overlapping elements. The proposed mesh enables the convergence of the mass, momentum and energy residuals criteria with acceptable orthogonal quality. Clustering techniques may be employed to calculate the mean view factor for grouped element grids instead of surface elements flow values.

[0096] A re-normalisation turbulence model ^^ − ^^ may be selected to model the air turbulence. The targeted cooling system in Figure 6 is a is Cold Air Distribution system, supplying air at relatively low temperatures for 10℃, to efficiently cool the zone with reduced energy consumption. The supply is simulated as an overhead inlet 0.3 meters in length × 0.3 meters in width with supply air velocity of 0.2 m / s, and the pressure outlet is 0.3 meters in length × 0.3 meters in width, with neutral pressure. A coupled scheme may be employed to solve the pressure-velocity coupling by brute force algorithms, for example by utilising SIMPLE, SIMPLEC, and PISO as segregated solving algorithms whereby the coupled algorithms offer robust single-phase implementation which can enable using larger mesh sizes for transient flows. Transient state time step size of the model is calculated based on the model minimum cell size of 11.7 mm in respect to the characteristic flow velocity of 0.2 m / s, starting at an initial condition at time ^^ = 0 and terminating at ^^ = 600. The initial phase may be discarded because of the slow start of zone conditions’ change. Thus, the dimensionality of the global state in the model is 25 × 108double precision arithmetic.

[0097] Thermal comfort PMV index’s empirical calculations are constituted by six factors: zone air temperature ^^^^, mean radiant temperature ^^^^, air velocity ^^^^and relative humidity ^^ ^^; besides clothing insulation ^^^^ ^^and human activity ^^. The values of airflow characteristics are provided by the Navier Stokes equations solutions of the fluent solver for each discretised cell, assuming dry bulb air temperature is the same as mean radiant temperature with a zone with no fenestration. The metabolic rate for human activity, garment insulation, and relative humidity are all fixed parameters for this example. Based on Equation (15), PMV index relationships is a polynomial quadratic equation, entailing the high dimensionality of the thermal comfort system and the importance of the proposed framework (described further with reference to Figure 8).

[0098] Set out below are User Defined Function (UDF) interpretation from a source file. The CFD model's output value for each time step is changed and the PMV values are saved. The PMV index representation may correspond for example, to the ISO 7730 thermal comfort categories (as described with reference to Figure 7). The outputof the CFD model may be presented as ASCII files (or similar) for each time step, creating the datasets for the autoencoder input layers. ^^ ^^ ^^ = [0.303 × exp(−0.036 ^^) + 0.0275] × { ^^ − ^^ − 3.05[5.733 − 0.007( ^^ − ^^) − ^^^^] − 0.42( ^^ − ^^ − 58.15) − 1.73 × 10−2× ^^(5.867 − ^^^^) − 0.0014 ^^(34 − ^^^^) −1.00 + 1.290 ^^ ^^ ^^ ^^ ^^ ≤ 0.078 ^^2. ^^ / ^^ ^^ =^^ ^^ ^^ ^^^^ ^^{ (18) 1.05 + 0.645 ^^^^ ^^^^ ^^ ^^ ^^^^ ^^> 0.078 ^^2. ^^ / ^^

[0099] In a further example, with reference to Figure 6, a cold air distribution system has the potential to maintain minimum thermal satisfaction according to ISO 773, while reducing the supply fan energy consumption. In this example, a CFD output is contributing as the neural network training, validation and testing sets. A CFD model may be applied where the geometry is simulating an indoor zone with dimensions of 3 meters in length × 3 meters in width × 2.6 meters in height. The CFD model mesh is unstructured with sensitivity analysis conducted to reach nearly 1,050,000 non overlapping elements. In order to reduce the computation requirements for training the model, the CFD mesh elements are reduced to a point cloud model with 600 points for 600 timesteps. The boundary conditions of the system are assigned to the inlet for velocity and temperature, heat load gains through the envelope and occupant 630 shown in Figure 6.

[0100] The boundary conditions are applied to a transient model solver using ANSYS Fluent R22.2. A User Defined Function (UDF) is implemented transiently through the model.

[0101] ^^ ^^ ^^ =[0.303 × ^^ ^^ ^^(−0.036 ^^)+ 0.028]×{^^ − ^^ − 3.05[5.733 − 0.007(^^ − ^^)− ^^ ^^)− 0.42(^^ − ^^ − 58.15)− 1.73 × 10−2× ^^(5.867 − ^^ ^^)− 0.0014 ^^(34 − ^^ ^^) − 3.96 × 10−8^^ ^^ ^^( ^^ ^^ ^^ + 273)4− ( ^^ ^^^^ ^^)}

[0102] Where M is the metabolic heat rate (W / m2) used as a constant for a standing human activity; W is the activity level (W / m2); ^^ ^^ ^^ is the clothes level temperature (℃) added as a constant for short sleeves shirt and trousers, Pa is the water vapor pressure extracted from the CFD simulation data; ^^ ^^ is the ambient air temperature is taken from the CFD solution, ^^ ^^ ^^ (Clo) is clothing thermal insulation; ^^ ^^ ^^ is the clothing ratio; ^^ ^^ is the radiant temperature (℃) considered equal to ^^ ^^ for simplification; ℎ ^^ is heat transfer convective coefficient (W m2 / ℃ ).

[0103] As for the control performance, the Koopman autoencoder approach is compared with other linearisation techniques performance with the application of LQR. The control input ^^^^is applied to the inlet velocity and temperature, with the desired state space ^^^^^^ ^^ ^^is defined where the PMV index ranges between -0.7 and 0.7 based on ISO 7730.

[0104] Figure 7 is a schematic diagram 700 illustrating PMV index variables and ISO7730 categories in accordance with an exemplary embodiment of the present invention. An occupants personal conditions, including metabolic rate and clothing insulation 705, is combined with environmental physical parameters of ambient temperature, radiant temperature, air velocity and relative humidity 710 through 715 and added to the convective transfer 720 and a clothing insulation ratio. The outcome of the calculation is based on PMV ISO 7730 categories 730, to obtain highest thermal comfort satisfaction, moderate satisfaction, minimum satisfaction or thermal discomfort. Network Architecture

[0105] Sequential feedforward networks are utilised with the present invention to exploit the consecutive order of the transient PMV-CFD model, passing the data input through all the modules in the same order defined through the linear transforming layers. In the autoencoder sequential network, each layer input is reshaped by a set of hidden layers represented as( ^^(1)( ^^))). The overall length of the hidden layers’ chain is equivalent to the depth of the model, while the node of each layer is equal to the width of it. Each linear hidden layer takes the form of ^^ ^^ + ^^ with ^^ representing the weight matrix and ^^ representing the bias matrix assigned to each layer. During autoencoder sequential training, the mapping of ^^( ^^) is driven to match ^^∗( ^^) with the nonlinearity accounted for using activation function, which enables theautoencoder to serve as a hierarchical features extractor. The parameters of the autoencoder are learned using gradient descent with respect to the loss function based on the different phases of our deep autoencoder. The loss function restricts the difference between the input and the reconstructed output, including a regularisation term to guarantee no redundancy in the features and reducing the impact of peripheral features on the network performance through weight regularization.

[0106] A challenge is presented for CFD mesh representation and mesh data, with standard neural networks designed to respond to pixelated equal spaced data from structured images. In addition, standard neural networks rely on cardinal directions of the image pixels whereas mesh can be given in any orientation. Several approaches tackled the treatment of mesh data, similar to pixelated image data such as: volumetric 3D boxes, implicit functions, mesh convolutional interpretation and point cloud data. In the CFD-PMV model which may be used with the present invention, point cloud approach is utilized with unordered X, Y, Z coordinates for its flexibility and reduction of the training dataset size. In the present invention, the size of the training is determined by multiplying the total number of point cloud samples by the quantity of time snapshots. Similar to images, the CFD-PMV mesh input for the network architecture has the first three dimensions represent the spatial dimensions, and the fourth dimension represents the airflow thermal comfort index conditions as shown in Figure 7. The point cloud sampling transforms the mesh from 1 million cell to 500 unequal spaced points as shown in Figure 8.

[0107] Figures 8a, 8b, 8c and 8d are point cloud diagrams illustrating PMV index point clouds for transient (a) c, (b) ^^^^= 50, (c) ^^^^= 200 and (d) ^^^^= 600 illustrating an indication of the nonlinear change of the PMV index distribution over 600 timesteps.

[0108] A machine learning framework may be utilised such that the autoencoder layers consist of linear sequential layers to address the transient state of the thermal conditions, with the depth of encoder network consisting of 5 layers and each layer contains parameters associated set of three-dimensional widths such as: 1024, 512, 64, 32, 8 and 8 in the six layers of the encoder. These layers are followed by an activation with the rectified liner unit ReLU: ^^( ^^) = ^^ ^^ ^^{0, ^^}, with experiments proving the faster performance of ReLU activation function in comparison with sigmoid. The weight matrix ^^ is randomly initialised from a Xavier uniform distribution in the rangeof [− ^^, ^^] for ^^ = 1 / √^^ where a is the input layer dimension. The bias vector is initialised to 0.

[0109] The auxiliary network is added to further reduce the PMV dimensionality, Koopman operator eigenvalues and eigenvectors as described with reference to Equations (1) to (9). The input to the auxiliary network is the encoded state spaces ^^̂, and it outputs the eigenvalues of ^^ as the complex conjugate pair ^^±= ^^ ± ^^ ^^, and + ^^̂22will be mapped to ^^ and ^^. This is achieved through two auxiliary networks denoting the one for the complex conjugate pair and another for real parts of the eigenvalues as mentioned in equation (9). For comparison, another autoencoder network is created with an alternative of auxiliary network applying SVD PCA method described with reference to Equations (10) to (13) provided directly as a known library (such as a python library in the context of the Pytorch framework for example), where the optimal encoder weights are going to be given by the pseudo-inverse of the decoder weights if the training data spans the full dimensionality. The multiplication of the^^and ^^^^^^may act as the hidden dynamic layer.

[0110] The performance of dimensionality reduction for both methods are further described with reference to Figures 9 to 13. With the objective to minimize the error, the explicit loss function of the autoencoder network follows the structure of the three phases as described with reference to Figure 5. The loss function has three mean error loss components: reconstruction loss ^^^^ ^^ ^^ ^^ ^^, parameterisation loss ^^^^and prediction loss ^^^^ ^^ ^^ ^^, followed by regularisation ^^2to the weight matrix. The total loss is calculated as in Equation (20):where ^^1, ^^2, ^^ ^^ ^^ ^^3are loss hyperparameters.

[0111] After defining the loss function and components, the training of the autoencoder is conducted for the training dataset, with an optimiser (ADAM for example) and an adaptive learning rate, storing an exponentially decaying average of past squares gradient descent. The hyperparameters of the network are tuned based on the validation loss performance for every epoch compared to the training losses as described with reference to Figures 9 to 13.

[0112] Advantageously, the Koopman operator continuous and discrete analysis combined with the autoencoder deep network reduces the nonlinearity of the PMV thermal comfort index. Figures 9 to 13 illustrate the operation and performance of the network across the training and validation datasets, whereby the results of the SVD network are shown as well as the eigenvalues of the SVD output for each training model are shown considering the discrete model stability through a unit circle (as will be described with reference to Figure 11). The performance of the SVD of the discrete model for the eigenvalues and eigenvectors is compared with the Koopman eigenfunction continuous model analysis avoiding a PMV index with high dimensionality.

[0113] In the present invention, autoencoder performance is optimised by utilising hyperparameter tuning. Figure 9 illustrates charts 905, 910, 915, 920, 925 and 930 for autoencoder hyperparameter tuning. The autoencoder is trained for a multiple layer configuration, batch sizes, optimiser learning rate, the trajectory timestep, and the losses hyperparameters, based on the lowest total validation loss as a metric shown in Figure 9. In this case, the results indicate that the encoder optimal structure is [501, 1024, 512, 64, 32, 8, 8], while the optimal layers’ structure for the decoder is [8, 16, 32, 64, 128, 501].

[0114] In addition, loss hyperparameters are optimised as ^^1= 0.0001, ^^2= 0.001, and ^^3= 0.001 as shown in charts 920, 925 and 930. Optimiser learning rate has been found to be the most significant parameter effecting the training speed - optimally reducing the validation error at 0.014 maximises the model performance. The discrete PMV index system is trained for 3 hours, for 200 epochs including early stopping where the lowest validation error is considered for each epoch. In this way, gradient vanishing is avoided with the least training epochs. This is shown in Figure 10 which is a chart 1000 illustrating network training versus validation performance through defined epochs in accordance with an exemplary embodiment of the present invention. As can be seen from Figure 10, the best training loss retrieved from the network is 0.021, matching a validation loss of 0.016.

[0115] Stability of the identified system ensures that it is under control and is a desirable characteristic for controllers’ design, specifically for the quadratic complex polynomial of PMV index. To investigate the stability of the latent space, the eigenvalues of the reduced dynamical PMV index determined how stable the system’s behaviour is. Mathematically, to attain a stable system, the eigenvalues should be close to a unit circle as shown in Figure 11.

[0116] Figure 11 is a chart 1100 illustrating eigenvalues λ_j visualisation for the real plane in accordance with an exemplary embodiment of the present invention. Since autoencoders can be difficult to train (when the eigenvalues of the hidden weight matrix deviate from the absolute value of 1), this may lead to the issue of vanishing and exploding gradients. Advantageously, overcoming this issue enables faster converging forecasts of the network. In the present invention, loss regularisation may be provided for nonconvex minimisation in loss calculations (which may be carried out in ‘Pytorch’ or the like). Figure 11 demonstrates that the reduced states ^^̂^^, ^^, ^^ = 1, … , ^^^^̂decomposed eigenvalues strictly converge within the unitary circle 1105. This indicates that the reduced states are dynamically stable, leading to improved predictions, following the least prediction error. The unitary matrices generated in this framework are parameterised in such a way that the gradient-based optimisation can be applied with the absolute values of eigenvalues <1 and preserve the hidden weight norms for powerful long-term dependencies for the sequential model of the present invention. Network Prediction Performance

[0117] Figure 12a is a PMV-CFD cloud point diagram 1200 illustrating the second time step in a testing trajectory dataset in accordance with an exemplary embodiment of the present invention. Figure 12b is a cloud point diagram 1205 showing Koopman autoencoder network cloud point predictions for an exact time step in accordance with an exemplary embodiment of the present invention. Figures 12a and 12b illustrate comparative visualisations of the PMV-CFD simulated results and reconstructed solutions showing the performance of the present invention. As can be seen in Figures 12a and 12b, the autoencoder model retains the important variation of the PMV index transient distribution and that that the general pattern of the PMV index distribution is maintained. In particular, the PMV index points are retained as are local characteristicsof occupancy in the room, inlet and outlet airflow distribution during the various the dimensionality reduction various stages.

[0118] Figure 13 is a chart 1300 illustrating PMV index prediction error relative to the CFD global state results with the Koopman parameterisation method and SVD method in accordance with an exemplary embodiment of the present invention. As can be seen from Figure 13, after the training process of the neural network, the PMV index based on CFD simulation testing data set is used to check the validity of the mode prediction. Each box 1305 represents the relative error of timesteps of the testing dataset trajectory and the predictions of the model. The approximate median error of the entire trajectory is less than 0.5%, with a minimum error of 0.01% and the highest maximum error of 2.5%. Outliers in the reconstructed data relative error to the CFD simulation data show that few points in the CFD simulated data possibly did not converge to a solution in the CFD ANSYS model. Advantageously, these results show high levels of accuracy for the PMV index reconstruction from Koopman basis and to provide reduced order modelling. Chart 1300 further illustrates a comparison between Koopman parameterisation and SVD methods. Chart 1300 represents the prediction performance of each method compared to the actual CFD readings. As can be seen, the Koopman parameterisation method achieves better prediction performance with less relative value percentages in each timestep in one trajectory with a minimum relative error of 0.001% and a maximum relative error of 2.23%(compared to SVD with a minimum of 0.01% and the highest maximum error of 2.5%). Further, Koopman parameterisation offers a linear model in the latent space, enabling linear feedback control application, while SVD only offers dimensionality reduction. Koopman parameterisation provides a superior method for identifying system nonlinearity when compared to the SVD method.

[0119] While the invention has been described in conjunction with a limited number of embodiments, it will be appreciated by those skilled in the art that many alternatives, modifications and variations in light of the foregoing description are possible. Accordingly, the present invention is intended to embrace all such alternatives, modifications and variations as may fall within the spirit and scope of the invention as disclosed.

[0120] Results: Autoencoder system identification with control input.

[0121] With the system and method of the present invention described above, it would be desirable to learn the entire dynamic behaviour of the system by only using the output of the autoencoder. As will be appreciated, controllability is essential to determine whether the system defined by the Koopman operator framework is capable of transitioning state ^^^^to state ^^^^+1. In an example, the encoder structure may transform the 1000-point cloud of CFD data based on the structure of [1000, 1024, 512, 64, 32, 8]. This structure may be fine-tuned by ‘Raytune’ as a hyperparameter to optimise the network performance. The decoder structure is defined as [8, 16, 32, 64, 128, 1000] and the defined controller input of the air supply velocity and temperature, the controller input dimension is defined as 1. Losses hyperparameters may also be fine-tuned and given the values of [0.0001, 0.001, 0.001, 4]. Other parameters, such as the batch size and learning rate of the ‘Adam’ optimiser may be set as 16 and 10 - 4, respectively. The output dimension 8 × 8 of the network acted as the dimension of the latent state space and the A matrix dimension with the values of A matrix composed by the dynamic weight. The B matrix dimension may be set as the dimension of the auxiliary network and the controller input as 8 × 1. Training for the network may be set up for different epochs ranging from 5 to 200 training epochs. It has been found that at the start of the training with 5 epochs, and applying the rank testing or the like, the results indicate that the system is uncontrollable with a high training error of 0.67. However, at the 200th epoch of training, it has been found that the system is controllable and observable by the full rank test (i.e. ^^ = [^^ ^^ ^^ ^^2^^ ⋯ ^^^^−1^^]) with a training error of 0.009 as shown in the chart 1400 of Figure 14 which plots epochs 1405 against total loss 1410 with plots of train_loss 1415 and val_loss 1420. As shown in chart 1400, there is a high validation error at the 110th epoch, which might be caused by the limited number of trajectories of the validation dataset. This implies the correlation between the training period of the autoencoder and the accuracy of the A and B matrices definition. While controllability and observability do not guarantee the reached state cannot be maintained by the controller they guarantee that the latent states given by the encoder and the auxiliary network represent a linear latent system

[0122] Results: Closed loop LQR performance of the linear latent system.

[0123] After ensuring the controllability of the Koopman framework system described above, it would be desirable to get the identified linear-latent system ^^ to the reference / desired state of PMV index neutrality and maintain that neutral thermalsensation throughout the control time horizon. For simplicity and in order to apply the LQR controller noted above, the Q cost matrix may be set up to perform as an identity matrix and the R controller regulator may be set up to be a vector of 1. These matrices may be deployed to calculate the initial K gain factor to assess the LQR performance and tune it accordingly. ^^0represents the initial condition of a CFD model with a supply inlet temperature and velocity of 10℃ and 0.2 m / s for a cold air distribution system. The previously trained encoder network with its weights and biases may be applied to the initial state ^^0and the reference / desired state ^^^^ ^^ ^^and ^^^^and ^^^^ ^^ ^^are calculated. The difference between the two states (based on ^^ = ^^^^^^ ^^ −(^^^^^^ ^^)(^^ + ^^^^^^ ^^)−1( ^^^^^^ ^^) + ^^) is multiplied to calculate the control signal u∗t. The u∗tmay be applied to the linear system of the latent space 1500 as shown in Figure 15 (and defined by ^^+= ^^ ^^ + ^^ ^^). The initial state is updated with the next ^^^^as well as the u∗t, for 500 timesteps. The performance of the control loop is shown in the chart 1600 of Figure 16 which is a chart of control timesteps 1605 against PMV Value 1610. The data 8 points on chart 1600 are the intrinsic coordinates given by the Koopman parametrisation. Each point represents the matrix of the lifted state ^^^^with the behaviour of each point as shown against the control horizon. The magnitude of these points are either ranging from very high or low PMV points indicating critical PMV values that require attention of the system. With the application of the closed-loop controller, all of the points are converging to the neutral PMV index with the calculated u∗tafter almost 250 timesteps, which demonstrates the effectiveness of the LQR controller on the latent space thermal performance. It will be appreciated that the period of reaching a neutral PMV index can be adjusted and reduced with tuning the Q and R matrices to reduce the controller effort.

[0124] Results: Closed-loop performance of the CFD model using the Koopman- based control law and a comparison with a RC-model-based control law.

[0125] The Koopman operator calculated controller input may be applied to several CFD models to test the performance of the LQR signal on the PMV index problem on the zone level rather than an 8-point control problem. This enables measuring the effectiveness of the derived controller input and the closeness to the desired PMV conditions; in this case, the desired outcome is to get the PMV CFD point cloud values as close as possible to 0. In the process of running the CFD simulations, ANSYS Fluentinlet boundary conditions starting at 10℃ and 0.2 m / s for the cold air distribution system may be used as the initial condition of the simulation. Different frequencies of applying the control input may be applied to examine whether the calculated Koopman operator can improve the PMV values. Figure 17A to 17D depicts four points as a sample of the CFD point cloud with the PMV values and timesteps application on the transient CFD model using the Koopman-model-based LQR controller K. Testing the control input for several frequencies shows that the longer the application of the controller input is, the higher the inlet velocity is and the closer the sample points are to the PMV values. Sample points given in Figure 18B show the different locations of each point. Points 100, 500, and 700 experience gradual improvement of the PMV values with an apparent decrease getting closer to the acceptable PMV category. It is evident that there is a strong correspondence exists between a decrease in the PMV values and the increase in the controller inlet velocity signal. This provides confidence that using Koopman operator with LQR controller is sensible for the improvement of the thermal comfort in the zone. However, the point 300 PMV value was steady and then increased due to the proximity to the occupant, impacted by the dissipation of the heat generated by the occupant in the CFD model. With the increase of the transient timesteps, the heat generated by the occupant would require an increase in the control signal velocity input and longer CFD simulations, hence demanding more computational resources. As the controller input is calculated by multiplying the ^^^^ ^^ ^^by the difference between the existing state and the desired state (as in ^^ = ^^^^^^ ^^ −(^^^^^^ ^^)(^^ + ^^^^^^ ^^)−1( ^^^^^^ ^^) + ^^) and can be scaled by this operator, the control inputs become quite interpretable. While these results are promising, further simulation horizons would be preferable to capture additional dynamics. Furthermore, these results are significant, as the controller input calculated by the Koopman framework indicated that it is capable of improving the PMV values without averaging the zone conditions. In the application of cold air systems, the Koopman operator control framework was capable of suppressing potential overcooling.

[0126] If we turn to the RC model as a system identification technique combined with PMV linearisation and LQR control signal calculations, the R value combined representing the enclosure of the zone is applied as 1.485 ^^^^. ^^ / ^^, the total thermal capacitance of 48.8 J / m2.K., and the outside air temperature is 30℃. A PMV-RC lumped LQR control strategy may then be applied in closed-loop to the RC lumpedmodel with the parameters ^^ = 100, ^^ = 1 and ^^ = 60 s, and PMV values are then measured over time – as shown in Figure 19 which is a chart 1900 of time 105 against PMV 1910 with plot 1915. It will be appreciated that an RC model including one averaged value for the zone will take a longer time to converge for the desired state.

Claims

The claims defining the invention are as follows:

1. A method for offline Koopman system identification, including the steps of: a) providing one or more generated snapshots ( ^^) representative of a thermal comfort quantification index, each being associated with a fixed range of boundary conditions for at least one of temperature, velocity ( ^^) or pressure (p); b) reshaping the one or more generated snapshots ( ^^) to generate one or more 3D tensors having variables (No. trajectory, ^^^^, ^^ ) and (No. trajectory, ^^^^, ^^) where: No. trajectory is the trajectory number; ^^^^is Timesteps for CFD horizon, sensor data time range or thermal images’ timesteps; ( ^^) are generated snapshots; and ( ^^) is a range of boundary conditions for temperature and velocity; c) splitting the generated 3D tensors into at least one training dataset, at least one validation dataset and at least one testing dataset; d) providing an encoder to reduce state space shape from ^^ to ^^^^= ^^^^( ^^^^); e) providing a dynamic auxiliary network having one or more layers to thereby calculate an A matrix and a B matrix, the A and B matrix being representative of a reduced order system definition; f) computing the auxiliary output of ^^̂^^+ ^^= ^^ ^^^^+ ^^ ^^ to … advance the reduced order system into a time step with the control of the boundary condition; g) providing a decoder such that the reduced order system may be reconstructed to the original state for a further time step from ^^̂^^+ ^^→ ^^̂^^+ ^^; and h) computing a loss function for reconstruction, linearisation and prediction.

2. The method of claim 1, further including training steps, wherein, over a predetermined number of epochs, for each of steps c) to h), the method further includes: i) calculating a loss function for the training dataset; j) calculating a loss function for the validation dataset; k) optimising the training process via an optimizer; and l) passing through stopping conditions of training loss vs validation loss.

3. The method of claim 1 or 2, wherein the method further includes the steps of: m) saving encoder hyperparameters A and B matrices; n) testing prediction accuracy against testing dataset for the autoencoder.

4. The method of claim 1, wherein the generated snapshots ( ^^) that are representative of a thermal comfort quantification index include Predictive Mean Vote - Computational Fluid Dynamic (PMV-CFD).

5. The method of claim 1, wherein the range of boundary conditions further includes relative humidity.

6. The method of claim 1, wherein at step b), reshaping includes transforming 3D meshes to a point cloud representation, thereby reducing computational resources.

7. The method of claim 1, wherein at step e), the dynamic auxiliary layer is created by way of:where: 1 ^̅^^^= ∑ ^^̂ , ^^ = 1, ... , ^^ , ^^^^ ^^̂^^and ^^ and ^^^^are orthogonal unitary colu ∑ ^^ ^^ ^^mns and^^= ^^ ^^ ^^ ^^( ^^^^, ^^)^^=1is a diagonal hierarchical matrix.

8. The method of claim 1, wherein at step e), the dynamic auxiliary layer is created by way of:where: ^^^^are eigenvalues represented by complex imaginary and real values ^^±= ^^ ± ^^ ^^ , denoting the oscillating and decaying nature of airflow and are parameterised by the cosine and sine matrix along with exponential representation.

9. The method of claim 1, wherein at step k), the optimizer includes one or more of an ADAM optimizer, Adadelta optimiser or SGD optimiser.

10. A method for generating an online feedback controller, including the steps of: a) embedding an encoder for ^^^^( ^^^^ ^^ ^^) → ^^^^ ^^ ^^, where: ^^^^, is a trained encoder network with optimal hyperparameters; b) applying a Linear-Quadratic Regulator (LQR) function on a predetermined set of Koopman operator matrices and a quadratic and regulator cost matrices to provide a K gain matrix; and c) setting the system environment to an initial condition ^^^^.

11. The method of claim 10, wherein at step a), the trained encoder network with optimal hyperparameters is the output of the offline Koopman system identification in claim 1.

12. The method of claim 10, wherein at step b), the predetermined set of Koopman operator matrices are the matrices A and B, derived from Claim 1.

13. The method of claim 10, wherein at step b), the K gain matrix is provided by the function ^^ ^^ ^^( ^^, ^^, ^^, ^^) → ^^ ^^ ^^ ^^.

14. The method of claim 10, wherein at step b), the K gain matrix is provided by an Model Predictive Control (MPC) function.

15. The method of claim 10, wherein for each t in timesteps for control ^^^^, the method further includes the steps of: d) embedding an encoder for the current state→ ^^; e) applying optimal feedback control ^^ ^^ ^^ ^^( ^^ − ^^^^ ^^ ^^); and f) adding disturbances update online.

Images