Method and system for solving large-scale ill-conditioned equation set for multi-energy network optimization

Abstract

The invention discloses a method and system for solving a large-scale ill-conditioned equation set for multi-energy network optimization, and the method comprises the steps: firstly building a multi-time-scale multi-target integrated energy system energy flow optimization model, and solving the multi-time-scale multi-target integrated energy system energy flow optimization model through a Newton-Raphson method; obtaining a large-scale ill-conditioned equation set; then carrying out block dimension reduction on a Jacobian matrix in the large-scale ill-conditioned equation set so as to decompose an iterative equation in the large-scale ill-conditioned equation set into iterative equations of all energy subsystems; and finally, decomposing the jacobian matrix subjected to block dimension reduction by using an incomplete LU decomposition method, and completing independent iterative solution of each energy subsystem by using a stable biorthogonal conjugate gradient method. The method is suitable for rapidly solving the large-scale ill-conditioned equation set when the multi-energy network is optimized, the model is reasonable, the algorithm is correct, economical efficiency and rapidity are met, and the method has engineering application value.

Description

Technical Field

[0001] The present invention relates to the field of integrated energy systems, and in particular to a method and system for solving a large-scale ill-conditioned equation group for multi-energy network optimization. Background Art

[0002] Integrated Energy System (IES) is an indispensable part of Energy Internet. Through multi-energy flow analysis and optimization, the system can realize the complementary utilization of various energy sources, coordinate energy conversion and use, improve energy utilization efficiency, promote regional economic and social development, and improve residents' living standards.

[0003] The multi-energy flow analysis method of integrated energy system is mainly used to study and optimize the interconnection and coordinated operation of different energy forms in the integrated energy system. Its core is to accurately describe the transmission, conversion and consumption process of gas, electricity, heat and other energy in the system by establishing a detailed mathematical model. By solving the model, the optimal operating state of the integrated energy system can be obtained for application in practical engineering.

[0004] Traditional single energy flow analysis methods usually only consider a certain form of energy. For example, the power system mainly focuses on power transmission and dispatching. However, with the diversification and integration of energy systems, the analysis method of a single energy form can no longer meet actual needs. Multi-energy flow analysis methods not only need to consider the independent operating characteristics of various energy sources, but also need to conduct in-depth research on the interactions and synergies between different energy forms. For example, in an integrated energy system that includes natural gas, electricity, and heat, natural gas can be used to generate electricity, electricity can drive heat pumps to generate heat energy, and waste heat can be transmitted to the user through the heat network for heating or cooling.

[0005] In addition to designing a multi-energy flow analysis method for an integrated energy system that includes gas / electricity / heat, domestic and foreign scholars are currently focusing on studying solutions with fast and high accuracy because the solution of the energy network equation of the integrated energy system has the characteristics of large-scale, nonlinear, and strong constraints. Some scholars have proposed optimizing the integrated energy system by decoupling the equipment layer and the low-level system for the problem of solving complex integrated energy systems of a single subject, and designed an individual-centered distributed parallel optimization method on this basis. However, regarding the problems of structural diversity, time scale differences, and complex coupling relationships between various parts in the energy flow network of the integrated energy system, there is still a lack of algorithms that can effectively improve the system's computational efficiency and accuracy.

[0006] When using the Newton-Raphson method model to perform multi-energy flow optimization calculations in an integrated energy system, a major problem is that there is a significant order of magnitude difference between network parameters and loads. This difference directly leads to poor condition numbers of the Jacobian matrix when applying the Newton method, which in turn makes its spectral distribution very dispersed, resulting in a pathological state of the matrix.

[0007] Specifically, a Jacobian matrix with a high condition number will cause the error propagation and amplification to become more serious in each iteration. This error accumulation effect will cause the Newton method to behave unstable during the solution process, causing the result after each iteration to deviate more and more from the true solution, which may eventually cause the algorithm to fail to converge to the correct optimal solution.

[0008] In addition, the dispersion of the Jacobian matrix spectrum will also have an adverse effect on the efficiency of solving linear systems. In each iteration of the Newton method, a linear system of equations needs to be solved. If the spectrum of the Jacobian matrix is ​​very dispersed, conventional linear solution methods such as direct inversion or Gaussian elimination will become very inefficient, further slowing down the overall calculation speed.

[0009] The existence of these numerical problems not only affects the convergence of Newton's method, but also has a negative impact on the speed and accuracy of the multi-energy flow optimization algorithm. The accuracy and speed of the optimization algorithm are important indicators to measure its performance, and these indicators directly determine the effectiveness and reliability of the integrated energy system in practical applications. If the optimization algorithm cannot quickly and accurately solve the multi-energy flow problem, the scheduling and management of the energy system will be severely restricted, and it will not be able to achieve its expected economic benefits and operational efficiency. Summary of the invention

[0010] In view of the defects in the above-mentioned background technology, the present invention proposes a method and system for solving large-scale ill-conditioned equations for multi-energy network optimization, which is suitable for multi-energy flow analysis and optimization of integrated energy systems. The method of the present invention can effectively solve the numerical problems caused by poor Jacobian matrix condition number, spectral distribution dispersion and ill-conditioned phenomena in Newton method calculation.

[0011] The present invention constructs a multi-time-scale multi-objective integrated energy system energy flow optimization model and uses the Newton-Raphson method to solve the model to obtain a large-scale ill-conditioned equation group. Considering the sparsity of the Jacobian matrix equation in the large-scale ill-conditioned equation group solved by the Newton-Raphson method, the block technology is used to reduce the dimension of the Jacobian matrix, so that the original iterative equation is decomposed into iterative equations of each energy subsystem, and then the incomplete LU decomposition method is used to decompose the reduced-dimensional matrix. On this basis, a subspace iterative processing method based on the Krylov-Arnoldi subspace algorithm is designed, and a stable biorthogonal conjugate gradient method is used to realize the separate iterative solution of each energy subsystem.

[0012] The objective of the present invention is achieved through the following technical solutions:

[0013] A method for solving a large-scale ill-conditioned system of equations for multi-energy network optimization, the method comprising the following steps:

[0014] Step 1: Establish a multi-time scale multi-objective integrated energy system energy flow optimization model, and use the Newton-Raphson method to solve the multi-time scale multi-objective integrated energy system energy flow optimization model to obtain a large-scale ill-conditioned equation group;

[0015] Step 2: performing block dimension reduction on the Jacobian matrix in the large-scale ill-conditioned equations, thereby decomposing the original iterative equations in the large-scale ill-conditioned equations into iterative equations of each energy subsystem;

[0016] Step 3: Use the incomplete LU decomposition method to decompose the Jacobian matrix after block dimensionality reduction, and use the stable biorthogonal conjugate gradient method to complete the separate iterative solution of each energy subsystem.

[0017] Furthermore, the step 1 specifically includes the following steps:

[0018] 1) Establishing a multi-time scale multi-objective integrated energy system energy flow optimization model, wherein the multi-time scale multi-objective integrated energy system energy flow optimization model includes a day-ahead multi-objective energy flow optimization model and an intra-day multi-objective energy flow optimization model;

[0019] Establish a day-ahead multi-objective energy flow optimization model.

[0020] In day-ahead scheduling, the objective function should take into account indicators such as daily operating costs, carbon emission costs and comprehensive energy efficiency, and formulate plans for the park's unit start-up and shutdown, energy storage charging and discharging, and electricity and gas purchases.

[0021] a) Daily operating costs include energy costs, operation and maintenance costs, and start-up and shutdown costs.

[0022] The energy cost expression is:

[0023]

[0024] Where: and are the electricity cost and natural gas cost under scenario s respectively; c PV Subsidize the price of photovoltaic power generation; is the electricity price in period t under scenario s, c GAS is the price of natural gas; and are the electricity purchase power and gas purchase power in period t under scenario s, is the air intake rate; T is 24, unit is h; Δt is the unit scheduling period.

[0025] The operation and maintenance cost expression is

[0026]

[0027] Where: i represents the type of equipment, including various energy conversion equipment and energy storage equipment, c M,i is the operation and maintenance cost per unit power output of device i; is the output power of device i.

[0028] The expression of start-stop cost is:

[0029]

[0030] Where: c SS,i represents the start-up and shutdown costs of equipment i, including gas turbines, gas boilers and electric hydrogen production systems; It is a 0-1 variable, indicating the operating status of device i, 1 for startup and 0 for shutdown.

[0031] b) Based on the carbon emission price, carbon emissions can be converted into carbon emission costs. The system carbon emission cost expression is:

[0032]

[0033] Where: F 2 is the carbon emission cost; s is the probability of scene s appearing; and are the equivalent carbon emission costs of purchasing electricity and gas respectively; is the carbon emission cost of the gas-fired units within IES.

[0034] c) Considering the differences between heterogeneous energy sources, Efficiency refers to comprehensive energy efficiency, and the comprehensive energy efficiency expression is:

[0035]

[0036] Where: F 3 For comprehensive energy efficiency; To consider the load power of energy quality; and They are respectively the energy storage charging and discharging power considering the energy quality; The purchased energy power taking into account the energy quality.

[0037] Establish a daily multi-objective energy flow optimization model.

[0038] Intraday real-time optimization is based on the day-ahead scheduling plan, taking into account indicators such as operating cost, carbon emission cost and adjustment cost.

[0039] The total cost expression of the scheduling period is:

[0040]

[0041] Where: F is the total cost of the daily dispatch period; t0 is the starting scheduling period; is the energy cost in period t; is the equipment operation and maintenance cost during period t; is the start-up and shutdown cost of the unit during period t; is the equivalent carbon emission cost of electricity purchased during period t; is the equivalent carbon emission cost of gas purchase during period t; is the carbon emission cost of the gas-fired units within IES during period t; Adjust costs for energy purchase plans; Adjust costs for unit output plans.

[0042] Adjustment costs include energy purchase plan adjustment costs and unit output plan adjustment costs.

[0043] The expression of energy purchase plan adjustment cost is:

[0044]

[0045] Where: and are the purchased power planned within the day and the day before, respectively; and are the gas purchase power planned within the day and the day before, respectively; and It is the adjustment coefficient for purchasing electricity and gas.

[0046] The unit output plan adjustment cost expression is:

[0047]

[0048] Where: k i is the output adjustment coefficient of equipment i; and are the output of equipment i in day-ahead and intraday scheduling, respectively.

[0049] 2) The Newton-Raphson method is used to solve the multi-time scale multi-objective integrated energy system energy flow optimization model to obtain a large-scale ill-conditioned equation set. The specific method is as follows:

[0050] Let the initial value of the state variable x be x (0) , then by Taylor expansion we have

[0051] F(x (0) +Δx (0) )≈F(x (0) )+J (0) Δx (0) =0 (9)

[0052] Where: J (0)For F(x) at x = x (0) The Jacobian matrix at this time; F(x) is the system deviation matrix, which can be obtained from the multi-time scale multi-objective integrated energy system energy flow optimization model; Δx (0) is the initial value correction matrix of the state variables.

[0053] Further available

[0054] F(x (0) )≈-J (0) Δx (0) (10)

[0055] Δx (0) ≈-[J (0) ] -1 F(x (0) ) (11)

[0056] Note Δx (0) =-[J (0) ] -1 F(x (0) ), the state variable is updated to

[0057] x (1) =x (0) +Δx (0) (12)

[0058] x (1) Use this as the new initial value and continue iterating to complete the calculation and solution.

[0059] Furthermore, the step 2 is specifically as follows:

[0060] Set the state variables of the power subsystem, natural gas subsystem, and thermal subsystem to x e =[U,θ] T , x g =[Π g ] T , x h =[Π h ,T s ] T , where U and θ are the amplitude and phase angle vector of the power subsystem voltage respectively; Π g is the vector of the square of the gas pressure at each node of the natural gas subsystem; Π h , T s are the pipeline flow and heating temperature vector of the thermal subsystem. Set the deviation matrix of the power subsystem, natural gas subsystem, and thermal subsystem as F e =[ΔP,ΔQ] T , F g =[ΔQ d ] T , F h=[ΔΦ,ΔT s ] T , where ΔP and ΔQ are the deviation vectors of active and reactive power in the power subsystem respectively; ΔQ d is the deviation vector of the natural gas flow rate at the natural gas subsystem node; ΔΦ and ΔT s are the deviation vectors of the node power and heating temperature of the thermal subsystem respectively. Then the Jacobian matrix after block division is

[0061]

[0062] Where: J ee , J hh , J gg The Jacobian matrices of the power subsystem, thermal subsystem, and natural gas subsystem are respectively; J eh , J eg is the block matrix of the power subsystem coupled with the thermal subsystem and the natural gas subsystem; J he , J hg is the block matrix of the thermal subsystem coupled with the power subsystem and the natural gas subsystem; J ge , J gh It is the block matrix of the natural gas system coupled with the power subsystem and the thermal subsystem respectively.

[0063] Note that the apparent power deviation vector in is the voltage phasor of the power subsystem; Y is the admittance matrix. Then we have

[0064]

[0065] Since ΔP and ΔQ are the real and imaginary parts of ΔS, respectively, the Jacobian matrix of the power subsystem is

[0066]

[0067] Taking partial derivatives of the natural gas subsystem, we can get the Jacobian matrix of the natural gas subsystem

[0068]

[0069]

[0070] Where: A gl is the node-branch association matrix of the natural gas subsystem branches.

[0071] Taking partial derivatives of the thermal subsystem, we can get the Jacobian matrix of the thermal subsystem

[0072]

[0073] Where: C his the specific heat capacity of the heat medium; T 0 It is the temperature of the heat medium just flowing out of the pipe, usually 40℃; A h1 b is the node-branch adjacency matrix of the thermal subsystem without the heat source node; b s , C s are respectively h The relevant coefficient matrix is ​​obtained from the network structure of the thermal subsystem.

[0074] The off-diagonal elements in the Jacobian matrix J are generated by the coupling between different energy subsystems, and the specific calculation method is as follows:

[0075] The regional integrated energy system is generally connected to the external large power grid. Considering that the CHP unit works in the mode of heat-to-electricity, the power fluctuation within the power subsystem is generally balanced by the large power grid, and will not affect the status of other subsystems in the integrated energy system. Therefore, there is J he =J ge =0.

[0076] In the natural gas subsystem, the natural gas network is connected to the external gas source. Similarly, when fluctuations occur within the natural gas subsystem, they will be balanced by the external gas source without affecting the status of other energy subsystems. Therefore, there is also J eg =J hg =0.

[0077] In the thermal subsystem, since the CHP operates in the mode of determining electricity by heat, the heat source of the thermal subsystem is the CHP unit. Therefore, when the internal thermal power in the thermal subsystem fluctuates, the thermal power output of the heat source will change, affecting the state of the power subsystem and the natural gas subsystem. eh and J gh is a non-zero matrix, and its expression is as follows:

[0078]

[0079]

[0080] Where: A hsource is the row of the heat source node in the thermal network node-branch association matrix; c CHP is the power distribution coefficient of the thermal power and electrical power output of the CHP unit; η CHP It is the gas-to-electricity energy conversion efficiency of the CHP unit.

[0081] In summary, the Jacobian matrix of the integrated energy system changes to:

[0082]

[0083] It can be obtained that the thermal system is relatively independent in the integrated energy system and can be solved separately to obtain the iterative step length Δx of this part. h After that, the corresponding power system and natural gas system can also be solved independently. That is, the calculation of the original iteration step Δx can be decomposed into:

[0084]

[0085]

[0086]

[0087] Where: Δx e is the iteration step of the power subsystem, Δx g is the iteration step size of the natural gas subsystem.

[0088] Furthermore, in step 3: in order to improve the convergence and stability of the iterative equations of each energy subsystem after Jacobian matrix block division, each iterative equation is preprocessed before solving. And the preprocessing method of the iterative equations of each energy subsystem is the same. Taking the iterative equations of the thermal subsystem as an example, the basic idea is to convert the linear equations into

[0089] J hh Δx h =-F h (25)

[0090] Transformed into an equivalent system of equations This makes the iterative algorithm converge faster.

[0091] A common method is to use a non-singular matrix K such that:

[0092] (J hh K -1 )(KΔx h )=-F h (26)

[0093] The incomplete LU decomposition method is used to decompose the Jacobian matrix after block dimension reduction, and the stable biorthogonal conjugate gradient method is used to complete the separate iterative solution of each energy subsystem. The specific method is:

[0094] First, the incomplete LU decomposition method is used to calculate the matrix K based on the Jacobian matrix of each energy subsystem, and then the stable biorthogonal conjugate gradient method is used to complete the iterative solution of the equivalent equations of each energy subsystem. For each energy subsystem, the Jacobian matrix decomposition and the solution process of the equivalent equations are the same. Therefore, taking the thermal subsystem as an example, the Jacobian matrix decomposition and the solution process of the equivalent equations are explained in detail. Specifically, the following steps are included:

[0095] 1) Using the incomplete LU decomposition method, the Jacobian matrix J of the thermal subsystem after block dimensionality reduction is hh Decompose and calculate the matrix K:

[0096] Considering that the Jacobian matrix of the thermal subsystem is generally a large asymmetric sparse matrix with dominant diagonals, the incomplete LU (ILU) decomposition method is a better choice. Therefore, the ILU decomposition method is used to decompose the Jacobian matrix of the thermal subsystem, and the decomposition method of the Jacobian matrices of the other subsystems is the same. In the ILU decomposition method, make J hh = Iteration matrix J in KR hh satisfy:

[0097] a) K = LU. Here L and U are not J hh It is not the LU decomposition of , but the LU decomposition of K.

[0098] b) Both L and U are sparse matrices

[0099] Its essence is an iteration. The following is an algorithm for constructing L, U, and R and a proof of the stability of the ILU method:

[0100] Assume J hh =(j ij ) n×n is an M-matrix, for any coordinate set We can always find a lower triangular matrix L with diagonal elements of 1 = (l ij ) and the upper triangular matrix U=(u ij ) and a residual matrix R = (r ij )satisfy:

[0101] l ij =0,if(i,j)∈P (27)

[0102] u ij =0,if(i,j)∈P (28)

[0103]

[0104] Make J hh =LU-R is a normal split, and L and U are unique. This theorem means that for any defined set of coordinates, there exists a normal split J hh = KR, so that the elements of the LU decomposition of the split matrix K on these coordinate sets are 0, thus ensuring the sparsity of L and U. The proof process is the construction process. Assuming that the size of the matrix is ​​n×n, it takes n steps to construct L, U, and R. Each step is similar to the Gaussian elimination method. Take the first two steps as an example:

[0105]

[0106]

[0107] Among them, L, R, and have special properties. R (1) is a matrix that only negates the coordinate elements that appear in P in the first row and the first column of the original matrix J hh , and the remaining elements are 0. The purpose is to make the corresponding coordinate elements in become 0; L (1) is a Gaussian elimination matrix, that is, a lower triangular matrix. Its purpose is to make the non-diagonal elements in the first column of all 0. Therefore, L(1) has no influence on the elements in the upper triangular part of . So, the upper triangular elements in are consistent with those in ; the elements of L(1) are also 0 at the coordinates specified in P.

[0108] After performing the above process for n steps, becomes an upper triangular matrix, and its elements at the specified positions in the set P are 0, obtaining the required U above. And from the iterative relationship, it can be known that:

[0109]

[0110] Note the characteristics of L and R. Since R (k) has non-zero elements only in the k-th row and the k-th column, while L (m) R (m) is to clear the elements below the main diagonal in the m-th column of the matrix R (m) . If the elements below the diagonal in the m-th column of the R (m) matrix are already 0, then L (m) = I, where I is the identity matrix. Therefore, when k < m, there is L (k) R (m) = R (m) . The above formula can be simplified to:

[0111]

[0112] Then L = (L (n-1) L (n-2) ...L (1) ) -1 , R = R (1) +…+R (n-1) , and there is J hh = K - R = LU - R.

[0113] 2) Use the stable bi-orthogonal conjugate gradient method to complete the iterative solution of the equivalent equations of the thermal subsystem.

[0114] Solve a system of linear equations:

[0115]

[0116] The general projection method is to use the Petrov-Galerkin condition to project from the m-dimensional affine subspace Δx 0 +K m (called search space) to find an approximate solution Δx m :

[0117]

[0118] Where: L m is another m-dimensional subspace (called constraint space); Δx 0 Represents the initial value of the solution iteration.

[0119] In the Krylov subspace method, K m It is the Krylov subspace:

[0120]

[0121] Where w is the residual

[0122] The Arnoldi algorithm is a method for processing asymmetric matrices projected onto the Krylov subspace. In the Arnoldi process, the recursive properties of the Krylov subspace are first noted:

[0123]

[0124] Assume K k-2 A unit orthogonal basis for K k-1 A unit orthogonal basis for So:

[0125]

[0126] so:

[0127]

[0128] After Schmidt orthogonalization, we can get:

[0129]

[0130] This can be rewritten as:

[0131]

[0132] Its matrix form is in is a (k+1)×k matrix, h i,j for The element in row i and column j, Only one subdiagonal line below the diagonal is non-zero, and the rest is equivalent to an upper triangular matrix.

[0133] Iterative Process of Krylov-Arnoldi Subspace Algorithm and Constraint Space L m It can be divided into orthogonal projection method, orthogonalization method and biorthogonalization method. It is asymmetric and highly sparse, so the iterative equation can be solved by using the biorthogonalization method based on stable biorthogonal conjugate gradients (BICGSTAB) To solve, the specific steps are:

[0134] a) Obtain the matrix K by following the calculation process in step 1).

[0135] b) Take Δx h The initial guess value is The allowable error is ε, which can be 0.001. make j=1.

[0136] c) Calculation If j-1 =0, then Exit the iterative process; otherwise go to step d).

[0137] d) If j = 1, let p j =r j-1 , otherwise let p j =r j-1 +β j-1 (p j-1 -ω j-1 λ j-1 ),

[0138] β j-1 =(ρ j-1 / ρ j-2 )(α j-1 / ω j-1 ).

[0139] e) By Solution calculate s j =r j-1 -α j λ j .

[0140] f) If ‖s j ‖≤ε, let Exit the iterative process; otherwise go to step g).

[0141] g) By Solution make ω j =(s j ,t j ) / (t j ,t j ),

[0142] h) If If the accuracy requirement is met, the iteration ends. Otherwise, let r j =s j -ω j t j , and add 1 to the value of j, and return to step c).

[0143] According to the calculated and Get Δx (j) , and Δx (j) Return to the step where the Newton-Raphson method solves the model to update the state variables and continue iterating until the solution is reached.

[0144] The present invention also provides a large-scale ill-conditioned equations solution system for multi-energy network optimization, which is used to execute the large-scale ill-conditioned equations solution method for multi-energy network optimization. The system comprises:

[0145] A modeling and preliminary solution module, for establishing a multi-time scale multi-objective integrated energy system energy flow optimization model, solving the multi-time scale multi-objective integrated energy system energy flow optimization model using the Newton-Raphson method, and obtaining a large-scale ill-conditioned equation group;

[0146] A block dimension reduction module is used to perform block dimension reduction on the Jacobian matrix in the large-scale ill-conditioned equation group, so as to decompose the iterative equations in the large-scale ill-conditioned equation group into iterative equations of each energy subsystem;

[0147] The decomposition and iterative solution module is used to decompose the Jacobian matrix after block dimensionality reduction using the incomplete LU decomposition method, and use the stable biorthogonal conjugate gradient method to complete the separate iterative solution of each energy subsystem.

[0148] The advantages of the present invention are:

[0149] The present invention takes the comprehensive energy system evaluation system such as environmental protection, energy efficiency, and economy as the optimization target, and proposes a method for constructing an energy flow optimization model of a multi-time scale and multi-objective comprehensive energy system. In view of the sparsity of the Jacobian matrix equation in the large-scale ill-conditioned equation group solved by the Newton-Raphson method, the Jacobian matrix is ​​pre-processed by dimensionality reduction using the block technology, so that the original iterative equation is decomposed into iterative equations of each energy subsystem. Subsequently, the matrix after block processing is decomposed using the incomplete LU decomposition method. On this basis, a subspace iterative processing method based on the Krylov-Arnoldi subspace algorithm is designed, and a stable biorthogonal conjugate gradient method is used to complete the separate iterative solution of each energy subsystem, thereby reducing the difficulty of energy flow calculation and shortening the energy flow calculation time. The present invention is suitable for quickly solving large-scale ill-conditioned equation groups during multi-energy network optimization. The model is reasonable, the algorithm is correct, and it meets the requirements of economy and speed, and has the value of engineering application. DETAILED DESCRIPTION

[0150] Example

[0151] The present invention is further described in detail below with reference to specific embodiments.

[0152] (1) Ill-conditioned matrix processing and convergence description

[0153] The convergence of an equation is related to its spectral distribution. If the eigenvalue of the iterative matrix is ​​large and relatively dispersed, it usually means that the equation does not converge. The spectral distribution of the coefficient matrix illustrates the effectiveness of the present invention in processing ill-conditioned matrices. Based on the IEEE standard system, the example was improved by adding a multi-energy coupling device for heating, cooling and electricity to the system for testing. Table 1 shows the maximum eigenvalue of the iterative matrix of different examples at a certain iteration in the energy flow calculation process before and after being processed by the method proposed in the present invention (including the block dimensionality reduction of the Jacobian matrix and the incomplete LU decomposition steps).

[0154] Table 1 Maximum eigenvalues ​​of the iterative matrix before and after processing

[0155]

[0156] It can be seen from Table 1 that after being processed by this method, the spectral distribution of the iterative matrix has been significantly improved. After optimization, the maximum eigenvalue of each example is less than 2, which is much smaller than the result before optimization.

[0157] (2) Comparison of the number of iterations

[0158] In order to verify the effectiveness of the method proposed in the present invention in reducing the number of iterations of the iterative method, multiple examples were selected for testing, among which the node 57 example to the node 300 example are from the IEEE standard, and the node 8971 example and the node 1648 example are obtained by splicing the balance nodes of multiple identical node 300 examples. On this basis, the above examples are improved by adding cold and hot electrical energy coupling equipment. Table 2 statistics the average number of iterations of different examples when solving the linear equation system when performing energy flow calculations, using the present invention, without optimization processing, and using Jacobi preprocessing, where the calculation accuracy is 0.001.

[0159] Table 2 The number of iterations when solving the linear equations in energy flow calculation using different processing methods

[0160]

[0161]

[0162] It can be seen that when performing energy flow calculation, the processing technology method of the present invention can significantly reduce the number of iterations and the calculation time. And the effect is better for large-scale systems, which can save more than 90% of the number of iterations. It is very suitable for online energy flow calculation of large-scale systems.

[0163] (3) Energy flow calculation acceleration analysis

[0164] In order to verify the effectiveness of the method proposed in the present invention in accelerating energy flow calculation, the present invention selected multiple examples for testing. The CPU used for calculation is Intel i7-4710MQ@2.5GHz, the memory is 16GB, the linear equations are solved by GPU general computing, the graphics card used is NVDIA GTX860M, and the examples used in the test are obtained by splicing multiple 300-node examples and adding cold and hot electrical energy coupling equipment.

[0165] Table 3 Comparison of energy flow calculation time consumption of different preprocessing methods

[0166]

[0167] It can be seen from Table 3 that the method of the present invention can greatly improve the energy flow calculation speed, and has great advantages over the commonly used Jacobi preprocessing method.

[0168] Another embodiment of the present invention further provides a system for solving a large-scale ill-conditioned equation system for multi-energy network optimization, the system being used to execute the method for solving a large-scale ill-conditioned equation system for multi-energy network optimization, the system comprising:

[0169] A modeling and preliminary solution module, for establishing a multi-time scale multi-objective integrated energy system energy flow optimization model, solving the multi-time scale multi-objective integrated energy system energy flow optimization model using the Newton-Raphson method, and obtaining a large-scale ill-conditioned equation group;

[0170] A block dimension reduction module is used to perform block dimension reduction on the Jacobian matrix in the large-scale ill-conditioned equation group, so as to decompose the iterative equations in the large-scale ill-conditioned equation group into iterative equations of each energy subsystem;

[0171] The decomposition and iterative solution module is used to decompose the Jacobian matrix after block dimensionality reduction using the incomplete LU decomposition method, and use the stable biorthogonal conjugate gradient method to complete the separate iterative solution of each energy subsystem.

Claims

1. A method for solving a large-scale ill-conditioned system of equations for multi-energy network optimization, characterized in that: The following steps are involved: Step 1: Establish a multi-time scale multi-objective integrated energy system energy flow optimization model, and use the Newton-Raphson method to solve the multi-time scale multi-objective integrated energy system energy flow optimization model to obtain a large-scale ill-conditioned equation group; Step 2: performing block dimension reduction on the Jacobian matrix in the large-scale ill-conditioned equations, thereby decomposing the iterative equations in the large-scale ill-conditioned equations into iterative equations of each energy subsystem; Step 3: Use the incomplete LU decomposition method to decompose the Jacobian matrix after block dimensionality reduction, and use the stable biorthogonal conjugate gradient method to complete the separate iterative solution of each energy subsystem.

2. The method for solving a large-scale ill-conditioned system of equations for multi-energy network optimization according to claim 1, characterized in that: The step 1 specifically includes the following steps: 1) Establishing a multi-time scale multi-objective integrated energy system energy flow optimization model, wherein the multi-time scale multi-objective integrated energy system energy flow optimization model includes a day-ahead multi-objective energy flow optimization model and an intra-day multi-objective energy flow optimization model; The objective function of the day-ahead multi-objective energy flow optimization model considers daily operating costs, carbon emission costs and comprehensive energy efficiency; a) Daily operating costs include energy costs, operation and maintenance costs, and start-up and shutdown costs; The expression of energy cost is: Where: and are the electricity cost and natural gas cost under scenario s respectively; c PV Subsidize the price of photovoltaic power generation; is the electricity price in period t under scenario s, c GAS is the price of natural gas; and are the electricity purchase power and gas purchase power in period t under scenario s, is the air intake rate; T is 24, unit is h; Δt is the unit scheduling period; The expression of operation and maintenance cost is: Where: i represents the type of equipment, including various energy conversion equipment and energy storage equipment, c M,i is the operation and maintenance cost per unit power output of device i; is the output power of device i; The expression of start-stop cost is: Where: c SS,i represents the start-up and shutdown costs of equipment i, including gas turbines, gas boilers and electric hydrogen production systems; It is a 0-1 variable, indicating the operating status of device i, 1 for startup and 0 for shutdown; b) The expression of carbon emission cost is: Where: F2 is the carbon emission cost; π s is the probability of scene s appearing; and are the equivalent carbon emission costs of purchasing electricity and gas respectively; is the carbon emission cost of the gas-fired units within IES; c) The expression of comprehensive energy efficiency is: Where: F3 is the comprehensive energy efficiency; To consider the load power of energy quality; and They are respectively the energy storage charging and discharging power considering the energy quality; To consider the purchased power of energy quality; The objective function of the intraday multi-objective energy flow optimization model considers the operating cost, carbon emission cost and adjustment cost, wherein the operating cost and carbon emission cost are the same as the daily operating cost and carbon emission cost calculation formula in the multi-objective energy flow optimization model a day ago; the objective function of the intraday multi-objective energy flow optimization model is specifically: Where: F is the total cost of the daily scheduling period; t0 is the starting scheduling period; is the energy cost in period t; is the equipment operation and maintenance cost during period t; is the start-up and shutdown cost of the unit during period t; is the equivalent carbon emission cost of electricity purchased during period t; is the equivalent carbon emission cost of gas purchase during period t; is the carbon emission cost of the gas-fired units within IES during period t; Adjust costs for energy purchase plans; Adjust costs for unit output plans; Adjustment costs include energy purchase plan adjustment costs and unit output plan adjustment costs; The expression of energy purchase plan adjustment cost is: Where: and are the purchased power planned within the day and the day before, respectively; and are the gas purchase power planned within the day and the day before, respectively; and It is the adjustment coefficient for electricity and gas purchase; The expression of the unit output plan adjustment cost is: Where: k i is the output adjustment coefficient of equipment i; and are the output of equipment i in day-ahead and intraday scheduling respectively; 2) The Newton-Raphson method is used to solve the multi-time scale multi-objective integrated energy system energy flow optimization model to obtain a large-scale ill-conditioned equation set. The specific method is as follows: Let the initial value of the state variable x be x (0) , then by Taylor expansion we have F(x (0) +Δx)≈F(x (0) )+J (0) Δx=0 (9) Where: J (0) For F(x) at x = x (0) The Jacobian matrix at this time; F(x) is the system deviation matrix, which is obtained by the energy flow optimization model of the multi-time scale multi-objective integrated energy system; Δx (0) is the initial value correction matrix of the state variables; Further F(x (0) )≈-J (0) Δx (10) Δx≈-[J (0) ] -1 F(x (0) ) (11) Note Δx (0) =-[J (0) ] -1 F(x (0) ), the state variable is updated to x (1) =x (0) +Δx (0) (12) x (1) Use this as the new initial value and continue iterating to complete the calculation and solution.

3. The method for solving a large-scale ill-conditioned system of equations for multi-energy network optimization according to claim 2, characterized in that: The step 2 is specifically as follows: Set the state variables of the power subsystem, natural gas subsystem, and thermal subsystem to x e =[U,θ] T , x g =[Π g ] T , x h =[Π h ,T s ] T , where U and θ are the amplitude and phase angle vector of the power subsystem voltage respectively; Π g is the vector of the square of the gas pressure at each node of the natural gas subsystem; Π h 、T s are the pipeline flow and heating temperature vector of the thermal subsystem respectively; the deviation matrices of the power subsystem, natural gas subsystem, and thermal subsystem are set as F e =[ΔP,ΔQ] T , F g =[ΔQ d ] T , F h =[ΔΦ,ΔT s ] T , where ΔP and ΔQ are the deviation vectors of active and reactive power in the power subsystem respectively; ΔQ d is the deviation vector of the natural gas flow rate at the natural gas subsystem node; ΔΦ and ΔT s are the deviation vectors of the node power and heating temperature of the thermal subsystem respectively; the Jacobian matrix after block division is: Where: J ee , J hh , J gg The Jacobian matrices of the power subsystem, thermal subsystem, and natural gas subsystem are respectively; J eh , J eg is the block matrix of the power subsystem coupled with the thermal subsystem and the natural gas subsystem; J he , J hg is the block matrix of the thermal subsystem coupled with the power subsystem and the natural gas subsystem; J ge , J gh The block matrices for coupling the natural gas subsystem with the power subsystem and the thermal subsystem respectively; Note that the apparent power deviation vector in is the voltage phasor of the power subsystem; Y is the admittance matrix; then ΔP and ΔQ are the real and imaginary parts of ΔS respectively, so the Jacobian matrix of the power subsystem is: Take partial derivatives of the natural gas subsystem and get the Jacobian matrix of the natural gas subsystem Where: A gl is the node-branch association matrix of the natural gas subsystem branch; Take partial derivatives of the thermal subsystem and get the Jacobian matrix of the thermal subsystem Where: C h is the specific heat capacity of the heat medium; T0 is the temperature of the heat medium just flowing out of the pipe; A h1 b is the node-branch adjacency matrix of the thermal subsystem without the heat source node; b s , C s are respectively h The relevant coefficient matrix is ​​obtained from the network structure of the thermal subsystem; The off-diagonal elements in the Jacobian matrix J are generated by the coupling between different energy subsystems, and the specific calculation method is as follows: Considering that the CHP unit works in the mode of determining electricity by heat, there is J he =J ge =0; In the natural gas subsystem, the natural gas network is connected to the external gas source, so there is J eg =J hg =0; In the thermal subsystem, since the CHP operates in the mode of heat-to-electricity, the heat source of the thermal subsystem is the CHP unit, so J eh and J gh is a non-zero matrix, and its expression is as follows: Where: A hsource is the row of the heat source node in the thermal network node-branch association matrix; c CHP is the power distribution coefficient of the thermal power and electrical power output of the CHP unit; η CHP is the gas-to-electricity energy conversion efficiency of the CHP unit; Then the Jacobian matrix of the integrated energy system changes to: In the integrated energy system, the thermal system is relatively independent and can be solved separately to obtain the iterative step length Δx of the thermal system. h After that, the corresponding power system and natural gas system can also be solved independently; that is, the calculation of the original iteration step Δx is decomposed into: Where: Δx e is the iteration step of the power subsystem, Δx g is the iteration step size of the natural gas subsystem.

4. The method for solving a large-scale ill-conditioned system of equations for multi-energy network optimization according to claim 3, characterized in that: In step 3, before decomposition, each iterative equation is first preprocessed, and the preprocessing method of all energy subsystem iterative equations is the same; among them, the preprocessing method of the thermal subsystem iterative equation is: The linear equations J hh Δx h =-F h (25) Transformed into an equivalent system of equations Make the iterative algorithm converge faster; Use a non-singular matrix K such that: (J hh K -1 )(KΔx h )=-F h (26)。 5. The method for solving a large-scale ill-conditioned system of equations for multi-energy network optimization according to claim 4, characterized in that: In step 3, the Jacobian matrix after block dimensionality reduction is decomposed by using the incomplete LU decomposition method, and the stable biorthogonal conjugate gradient method is used to complete the separate iterative solution of each energy subsystem, wherein the decomposition and solution method of the Jacobian matrix of each energy subsystem is the same; for the Jacobian matrix J of the thermal subsystem hh , specifically including the following steps: 1) Use the incomplete LU decomposition method to decompose the Jacobian matrix of the thermal subsystem after block dimensionality reduction: Jacobian matrix J of thermal subsystem hh The decomposition steps are as follows: During the decomposition process, J hh = Iteration matrix J in KR hh satisfy: a) K = LU, and L and U are not J hh LU decomposition of , but LU decomposition of K; b) Both L and U are sparse matrices; Assume J hh =(j ij ) n×n is an M-matrix, for any set of coordinates There exists a lower triangular matrix L with all diagonal elements being 1 = (l ij ) and the upper triangular matrix U=(u ij ) and a residual matrix R = (r ij )satisfy: l ij =0,if(i,j)∈P (27) u ij =0,if(i,j)∈P (28) Make J hh =LU-R is a regular split, and L and U are unique; Assuming the size of the matrix is ​​n×n, constructing L, U, and R requires n steps, each of which is similar to Gaussian elimination; the first two steps are: Among them, R (1) is to only transform the original matrix J hh The coordinate elements in the first row and the first column of P are reversed, and the remaining elements are 0. The corresponding coordinate elements in become 0; L (1) is the Gaussian elimination matrix, i.e., the lower triangular matrix, the purpose is to make The off-diagonal elements in the first column of are all 0; The upper triangular elements in The elements of L(1) are also 0 at the coordinates specified in P. After the above process is carried out for n steps, That is, it is an upper triangular matrix, and its specified position element in the set P is 0, and the required U is obtained; and from the iteration relationship, it can be seen that: When k < m, there is L (k) R (m) = R (m) , Equation (34) is simplified to: Then L=(L (n-1) L (n-2) ...L (1) ) -1 , And there is J hh =KR=LU-R; 2) The stable biorthogonal conjugate gradient method is used to iteratively solve the equivalent equations of the thermal subsystem: Solve a system of linear equations: The specific solution steps are as follows: a) Obtain the matrix K according to the calculation process of step 1); b) Take Δx h The initial guess value is The allowable error is ε, and the calculation make j = 1; c) Calculation If j-1 =0, then Exit the iterative process; otherwise go to step d); d) If j = 1, let p j =r j-1 , otherwise let p j =r j-1 +β j-1 (p j-1 -ω j-1 λ j-1 ), b j-1 =(ρ j-1 / r j-2 )(a j-1 / h j-1 ); e) By Solution calculate s j =r j-1 -α j λ j ; f) If ‖s j ‖≤ε, let Exit the iterative process; otherwise go to step g); g) By Solution make ω j =(s j ,t j ) / (t j ,t j ), h) If If the accuracy requirement is met, the iteration ends. Otherwise, let r j =s j -ω j t j , and add 1 to the value of j, and return to step c); The calculated and Get Δx (j) , and Δx (j) Return to the step where the Newton-Raphson method solves the model to update the state variables and continue iterating until the solution is reached.

6. A system for solving large-scale ill-conditioned equations for multi-energy network optimization, used to execute the method according to any one of claims 1 to 5, characterized in that: include: A modeling and preliminary solution module, for establishing a multi-time scale multi-objective integrated energy system energy flow optimization model, solving the multi-time scale multi-objective integrated energy system energy flow optimization model using the Newton-Raphson method, and obtaining a large-scale ill-conditioned equation group; A block dimension reduction module is used to perform block dimension reduction on the Jacobian matrix in the large-scale ill-conditioned equation group, so as to decompose the iterative equations in the large-scale ill-conditioned equation group into iterative equations of each energy subsystem; The decomposition and iterative solution module is used to decompose the Jacobian matrix after block dimensionality reduction using the incomplete LU decomposition method, and use the stable biorthogonal conjugate gradient method to complete the separate iterative solution of each energy subsystem.

Images