Automated Design Method for Digital Metasurface Cells Based on Characteristic Mode Theory

Abstract

An automated design method for digital metasurface cells based on the characteristic mode theory includes the following steps: Step 1: Preset the basic structure of the cell; Step 2: Preset the optimization region for genetic iteration based on the basic structure obtained in Step 1 and perform grid processing on the metal patch layer; Step 3: Set the population size POP and the number of iterations GEN of the genetic algorithm for each grid individual in the metal patch layer after grid processing; Step 4: Set the objective function of the algorithm for the unit structure corresponding to each grid individual in the population obtained in Step 3; Step 5: Solve the characteristic modes of the grid structures representing the metal patch layer structure in the population and calculate the objective function; Step 6: Generate offspring population individuals through genetic iteration of the objective function obtained in Step 5; Step 7: Obtain the optimal individual, and the automated design process ends. The present invention presets a clear objective function according to the desired physical characteristics of the cell and completes the automated design of the unit structure with the desired physical electromagnetic characteristics.

Description

Technical Field

[0001] The present invention belongs to the field of electromagnetic technology, and particularly relates to an automated design method for digital metasurface units based on the characteristic mode theory. Background Art

[0002] Currently, passive units cannot achieve the function of flexibly and real-time controlling the beam, so novel bit-type units are proposed. Currently, the design methods for novel electronically controlled units are as follows: ① Change the bias voltage applied across the varactor diode, and study and analyze the variation of the resonance characteristics of the metasurface unit under different voltages to achieve continuously adjustable reflection wave phase or a phase difference of 180°; ② Change the patch structure parameters loaded with PIN diodes, study the influence of each structure parameter on the reflection wave phase, and obtain the desired reflection wave phase control performance through structural parameter adjustment.

[0003] Although the existing novel electronically controlled bit-type units solve the problem of the flexibility of array control and can achieve multifunctional applications of the metasurface by real-time switching of the array coding through a single-chip microcomputer or FPGA, the existing unit design methods are all based on simulating the reflection coefficient of the unit, and the designer completes the process of simulating parameter scanning and optimizing the structural dimension parameters. For the designer, this process requires multiple unnecessary and time-consuming repetitive tasks such as comparison and adjustment.

[0004] When designing a reflective phase modulation unit loaded with diode elements, each time the structural parameters are adjusted, the unit characteristics under two equivalent circuit forms of diode conduction and truncation need to be considered separately. Therefore, during the structural design process, it is often necessary to repeatedly compare the unit effects and adjust the structural parameters, making the unit characteristic evaluation process more complex, resulting in the problem of over-reliance on the design experience of antenna designers and the repetitive parameter adjustment and optimization process. Summary of the Invention

[0005] In order to overcome the deficiencies of the above-mentioned existing technologies, the purpose of the present invention is to provide an automated design method for digital metasurface units based on the characteristic mode theory, aiming to enable the designer to only preset a clear objective function according to the desired physical characteristics of the unit and use the means of genetic algorithm iterative optimization to complete the automated design of the unit structure with the desired physical electromagnetic characteristics, without relying on manual adjustment, but directly completing the structural design with the desired characteristics.

[0006] In order to achieve the above purpose, the technical solution adopted by the present invention is:

[0007] An automated design method for digital metasurface units based on the characteristic mode theory, comprising the following steps;

[0008] Step 1: Preset the basic structure of the digital metasurface unit;

[0009] Step 2: Use the genetic algorithm to preset the optimization region for genetic iteration according to the basic structure obtained in Step 1 and perform grid processing on the metal patch layer;

[0010] Step 3: Set the population size POP and the number of iterations GEN of the genetic algorithm for the grid-processed metal patch layer;

[0011] Step 4: Set the objective function of the algorithm for each unit structure represented by the grid individuals in the population obtained in Step 3;

[0012] Step 5: Solve the eigenmodes of the grid structures representing the metal patch layer structure in the population and calculate the objective function;

[0013] Step 6: Generate offspring population individuals through genetic iteration of the objective function obtained in Step 5;

[0014] Step 7: Obtain the optimal individual, and the structural automated design process ends.

[0015] The basic structure in Step 1 is a reflective digital metasurface unit structure loaded with PIN diodes, which is a basic structure unit composed of a first metal layer 1, a first dielectric layer 2, and a second metal layer 3. According to the center frequency point f of the expected operating frequency band of the metasurface unit structure, the corresponding wavelength is λ, and the range of the unit structure period P can be preset as and the PIN diode model is selected. In the present invention, the first metal layer 1 is a metal patch layer, the second metal layer 3 is a full-coverage metal layer under the unit structure period, and the first dielectric layer 2 is a dielectric substrate material.

[0016] Furthermore, the step 2 takes into account the PIN diode model selected in the step 1, presets the position for adding the diode, and sets the diode introduction position to the center of the structure by means of a slit loading method in the center of the first metal layer 1; then presets the desired design form of the basic structure of the unit, adopts a multi-resonant slit structure to ensure the working bandwidth of the unit, and adopts a rectangular metal sheet containing a "艹"-shaped slit as the initial structural form of the first metal layer 1. By gridding the metal patch layer, its initial structure is cut into the form of countless identical small grids of Nx rows and Ny columns, wherein the "艹"-shaped slit is the optimized area in the first metal layer 1 where the grid needs to be randomly retained and deleted, and adopts a The grid is divided into small square grids, and the grid side length ml should be less than one twentieth of the wavelength λ at the center frequency point f. The grid size is as small as possible, which will make the genetic optimization algorithm contain more gene bits. In the present invention, the first metal layer 1 is gridded to facilitate the subsequent operation of step three. The small grids in the optimization area preset in step two are randomly coded and assigned 0 or 1 to indicate whether the small grid is retained. The assignment of 0 indicates that the small grid is deleted, and the assignment of 1 indicates that the small grid is retained. The grids outside the preset optimization area of ​​the first metal layer 1 are always coded and assigned 1, that is, the grids are always retained, so as to achieve the purpose of completing the structural optimization only in the predetermined optimization area.

[0017] Furthermore, in step three, the number of individuals in the initial population is set to POP, and an initial population consisting of POP unit structure grid individuals with different structures is obtained by generating POP random binary arrays; the number of iterations is set to GEN, and the more POPs in the initial population means the more individuals that can be selected from the parents, and the larger the number of iterations GEN means the closer the final optimization accuracy is to the target.

[0018] Furthermore, the research object of the objective function to be calculated in step 4 is the unit structure form corresponding to the representation of each grid individual in the population obtained after steps 1, 2, and 3, and the characteristics of the objective function are related to the electromagnetic physical properties (characteristic mode phase difference) of the final desired digital unit structure. The objective function set in step 4 is to complete the objective function calculation under the structural form represented by each grid individual in the population in step 5. The specific settings of the two objective functions of genetic optimization in step 4 are as follows:

[0019] Objective function function 1: The polarization mode of the incident wave is specified as Y polarization or X polarization, and the characteristic main mode polarization corresponding to the metasurface unit structure grid in the on-off state is required to be the same as the incident wave. The target polarization angle polar should tend to 90° or -90°, as shown in formula (1):

[0020] function 1=|polar|-90 (1)

[0021] Objective function 2: The digital phase difference of the characteristic modes is determined by the superposition of the phase difference of the expansion coefficients of the main mode modes in two states of the electronic control unit under the same incident wave, the difference in characteristic angles, and the difference in polarization angles. Therefore, the desired characteristic mode digital phase condition for designing an electronic control digital unit with a reflection field phase difference of Δgoal phase is as follows:

[0022] function 2 = ΔMWC + ΔCA + Δpolar - Δgoal phase (2)

[0023] As described in Equation (2), where ΔMWC is the phase difference of the expansion coefficients of the main mode modes in the on-off two states, ΔCA is the difference in the main mode characteristic angles between the two, Δpolar is the difference in the main mode polarization angles between the two, and Δgoal phase is the desired phase difference of the design goal. The optimization process of minimizing the optimization goal makes function 2 tend to 0. With iterative optimization, the characteristic mode digital phase difference characteristic of the target digital unit will approach the desired phase difference of the design goal.

[0024] Furthermore, Step 5 is based on the desired optimization objective function represented by the relevant characteristic mode physical parameters (main mode expansion coefficient, characteristic angle, polarization angle) that can reflect the characteristic field phase in Step 4. Combining with the NSGA-ii genetic optimization algorithm, the script is written using the MATLAB language, and the model establishment, parameter solution, and data result recording are completed by calling the electromagnetic simulation software CST Studio Suite 2020;

[0025] 5-1) Model establishment

[0026] Generate POP random binary arrays according to Step 3. Each array is in the form of a discrete grid of individual structures represented by random binary information. The POP grid individuals form the initial population. Sequentially read the binary information of the POP grid individuals and call the electromagnetic software CST Studio Suite 2020 to complete the construction of the structural model at the corresponding positions of the grid;

[0027] First, call the CST software and create a new project file. Considering that each grid individual in the initial population is a random binary row vector array, this row vector array is converted into an array matrix with Nx rows (total number of grid rows) and Ny columns (total number of grid columns). Then, read the array matrix row by row, and determine whether the grids within the optimization region are retained according to the binary information (grids outside the optimization region are all retained as described in step two). If retained, construct a grid structure model at the corresponding position until the construction of the overall upper-layer metal grid structure corresponding to the entire array matrix is completed. Then, add a new medium and set the dielectric constant according to the parameters of the intermediate dielectric layer board, and complete the construction of the intermediate dielectric layer structure according to the unit period P and the dielectric thickness h. Finally, complete the construction of the fully covered metal floor structure on the lower surface of the dielectric layer according to the unit period size. Thus, the construction of the metasurface unit structure based on one grid individual within the population is completed;

[0028] 5-2) Parameter Solving

[0029] For the metasurface units represented by each grid individual, after completing the copy and translation operations according to the unit period P and arranging them into a 4×4 finite periodic structure, set the field monitors at the frequency f: far field farfield, magnetic field h-field. Then, use the high-frequency integral equation solver as the engineering solver (this engineering solver supports eigenmode solving and sets the number of eigenmode patterns ModeNum for the solution analysis). Next, complete the setting of the Y-polarized plane wave excitation source incident perpendicularly to the unit structure along the -Z axis, and set the electric field vector ’EVector’ to (0,1,0). If it is desired that the plane wave excitation source is an X-polarized wave incident perpendicularly, just change the electric field vector ’EVector’ to (1,0,0). Additionally, considering that the unit structure is a reflective metasurface unit and the lower surface of the unit is a fully covered metal floor, set the engineering model boundary condition Zmin to ’electric’. Thus, the solution setting of the project is completed, and the engineering eigenmode solution is started;

[0030] 5-3) Data Result Recording

[0031] When the project completes the solution calculation, it is necessary to record the eigenmode parameter data used for the objective function calculation, including the amplitude and phase of the mode expansion coefficients, eigenangles, and eigenfields. This process is achieved by using MATLAB to call CST to read the project results and record the data. Record the far field, eigenangles, and the amplitude and phase data of the mode expansion coefficients of the eigenmodes successively according to the eigenmode ordinal number, and export the data results in ASCII format according to the project path and name them successively according to the mode ordinal number, and store them in each parameter folder;

[0032] So far, the establishment of the unit structure model, the solution of the characteristic mode parameters, and the recording of the data results have been completed, and the solution of the characteristic modes of the structural individuals in the population in the fifth step has been completed.

[0033] In order to achieve the calculation objective function of the fifth step, it is necessary to judge the dominant mode of the characteristic mode and complete the extraction of the relevant data of the objective function. The dominant mode of the characteristic mode has two characteristics: the mode expansion coefficient is relatively large; the polarization mode of the characteristic mode is the same as the polarization mode of the excitation source. Based on the above two characteristics as the judgment basis, the specific process is as follows:

[0034] First, read the mode expansion coefficient MWC saved in the previous step, find the characteristic mode ordinal number with the largest amplitude of the mode expansion coefficient, and read and record the phase of the mode expansion coefficient of this mode ordinal number as MWC_phase; then read the X-component amplitude |Ex| and Y-component amplitude |Ey| in the far-field component data of this mode ordinal number. If it is excited by a Y-polarized wave, the mode polarization angle is recorded as: polar = tan(|Ey| / |Ex|), and when it is excited by an X-polarized wave, the mode polarization angle is recorded as: polar = tan(|Ex| / |Ey|); then read the characteristic angle of this mode ordinal number and record it as CA. For each grid individual in the population in the third step, there are two structural forms under the on-off state of the diode. The difference between the two is only the presence or absence of the grid structure at the diode position in the engineering structure model. Therefore, the two can respectively use the above steps to complete the solution and recording of the characteristic mode parameters in the two states of the on-off of the diode. Therefore, the polarization angle polar for the objective function 1, the phase difference ΔMWC of the mode expansion coefficient, the difference ΔCA of the characteristic angle, and the difference Δpolar of the polarization angle for the objective function 2 can be calculated, so as to complete the calculation of the two objective functions, and record the objective function results in the binary gene bit information of the grid individual, which is convenient for the subsequent step six to select the optimal offspring population according to the objective function value.

[0035] Furthermore, in the sixth step, according to the objective function values recorded in the individual gene bit information, non-dominated sorting and crowding degree calculation are completed. The specific process is as follows:

[0036] 6-1) Non-dominated sorting

[0037] According to the objective function values, the individuals in the population are stratified by non-dominated sorting. The lower the non-dominated level number to which an individual belongs, the lower the priority and the higher the dominant position. The specific process of non-dominated sorting is as follows:

[0038] (1) The population size is N, and the objective function values of the individuals in the population are x(1), x(2)... x(N);

[0039] (2) Let i = 1;

[0040] (3) For all j = 1, 2,..., N and j ≠ i, the dominance relationship between individual i and individual j is determined by comparing the objective function x(i) of the i-th individual and the objective function x(j) of the j-th individual. If there is no x(j) that is superior to x(i), then x(i) is recorded as a non-dominated individual;

[0041] (4) Let i = i + 1, and repeat to (2) until all non-dominated individuals are found.

[0042] The set of non-dominated individuals obtained through the above processes (1) to (4) serves as the first non-dominated layer (priority 0) of the population. Then, the other individuals in the population except the first non-dominated layer are repeatedly processed through (1) to (4) to obtain the second non-dominated layer (priority 1), and so on until the dominance relationship of the entire population is hierarchically sorted;

[0043] 6-2) Crowding degree calculation

[0044] Calculating the crowding degree is to estimate the crowding degree of other solutions around a specific individual solution in the population. The crowding degree is the average distance between adjacent individual solutions calculated according to each objective function. The specific process of calculating the crowding degree is as follows:

[0045] (1) Sort the population individuals in ascending order according to the magnitude of each objective function value;

[0046] (2) For each objective function, perform non-dominated sorting on the population. Take the individual solutions corresponding to the maximum objective function value fmax and the minimum objective function value fmin as the boundary solutions, and set the crowding degree of the boundary individual solutions to an infinite distance. The crowding degree of other intermediate individual solutions is the normalized absolute difference of the function values of the adjacent individual solutions on both sides. Denote the m-th objective function value of the j-th individual solution as f (m) (j). In the case where this individual solution does not belong to the boundary solution, its crowding degree value is expressed as (f (m) (j + 1) - f (m) (j - 1)) / (f (m) max - f (m) min);

[0047] After non-dominated sorting and crowding degree calculation, each individual in the population has these two attributes: non-dominated layer (i.e., priority) and crowding degree value, and these are used as the basis for comparing the quality of individuals in the population;

[0048] If the non-dominated layer to which the i-th individual x(i) belongs is superior to the non-dominated layer to which the j-th individual x(j) belongs, then individual x(i) is superior to x(j); if x(i) and x(j) belong to the same non-dominated layer, then it is necessary to compare their crowding degrees. If the crowding degree f (m) (i) of the m-th objective function of x(i) is greater than the crowding degree f (m)(j), then x(i) is superior to x(j). Therefore, according to the non-dominated sorting sequence size and crowding degree of individuals, individuals are preferentially selected to form a new population for subsequent iterations. Repeat step five for the individuals in the new population until the objective function is optimal or until the maximum number of iterations is reached and stop.

[0049] Furthermore, the optimal individual is obtained in step seven, and the process of the genetic optimization algorithm ends. The on-off two states corresponding to the obtained optimal individual are the digital units that can be automatically designed and used as reflective phase shifters. In step two above, the addition positions of the PIN diodes are specified. Just introduce the equivalent circuit of the diode being on or off here directly. Thus, based on the characteristic mode phase condition, the automated design of the encoded metasurface unit loaded with diodes is completed. The electromagnetic characteristics of the digital unit formed after the optimal structure designed automatically is in the form of the equivalent circuit of the loaded diode can meet the characteristic mode digital phase at the central operating frequency f and conform to the design conditions.

[0050] Advantages of the present invention:

[0051] The present invention first derives the digital phase condition based on the characteristic mode theory as the technical principle for designing metasurface units loaded with active electronically controlled components. Based on the physical characteristic parameters of the metasurface unit structure reflected by the characteristic mode theory, the physical characteristics of the optimization target are made clearer and are reflected as the convergence target characterized by clear electromagnetic characteristics in the optimization algorithm. This automated optimization design combined with the theoretical parameter concept makes the final design result more meet the design expectation target and improves the quality of the optimization design.

[0052] The technical solution used in the present invention solves the problem that the design process of metasurface units relies too much on design experience and repeatedly adjusts and optimizes structural parameters, and improves the design efficiency. Based on the structural form, phase, frequency band and other conditions set by the designer according to the expected target function, the structural design is directly completed. The designer only needs to adjust the preset optimization region in step two of the technical solution of the present invention, the number of population individuals and the number of iterations that determine the optimization process in step three, and the preset objective function in step four. During other processes, there is no need for manual structural adjustment and optimization, and the automated design of the encoded metasurface unit structure with certain expected physical characteristics can be completed based on the expected basic form of the unit. This greatly reduces the workload of the designer, shortens the design process, and improves the design efficiency. Brief Description of the Drawings

[0053] Figure 1 is the flowchart of the iterative optimization of the genetic algorithm adopted by the present invention.

[0054] Figure 2 is the three-dimensional layer schematic diagram of the preset basic unit structure in step one of Embodiment 1 of the present invention.

[0055] Figure 3 It is a schematic diagram of the optimized area of the preset genetic iteration in Step 2 of Embodiment 1 of the present invention.

[0056] Figure 4 It is the on-state structure of 1-bit cells corresponding to the optimal individual grid finally obtained in Step 7 of Embodiment 1 of the present invention.

[0057] Figure 5 It is the cut-off state structure of 1-bit cells corresponding to the optimal individual grid finally obtained in Step 7 of Embodiment 1 of the present invention. Detailed implementation manners

[0058] The present invention will be further described in detail below with reference to the accompanying drawings.

[0059] The present invention aims to solve the technical problem of realizing the automated design of digitalized encodable metasurface units with inversion characteristics from the perspective of the bit-type digitalized phase conditions characterized by the physical parameters of the eigenmode theory, based on the physical characteristics of the metasurface unit structure reflected by the eigenmode theory. In view of the scattering problem of the structure body, the relevant parameters of the eigenmode theory clearly characterize physical concept information such as modal current, modal far field, main excitation mode expansion coefficient, eigenangle, etc. Therefore, from the perspective of the eigenfield phase, the digitalized phase condition of the eigenmode theory is derived, and this digitalized phase condition is combined with the genetic optimization mechanism to realize an algorithm with this phase condition as the optimization design objective function, so as to achieve the automated design of digitalized units. This optimization design process is completely based on the physical characteristics of the structure, so the finally achieved electromagnetic characteristics are more in line with the expected design goals, and the purpose of simplifying the structure design process and shortening the optimization process is achieved.

[0060] The purpose of the present invention is to elaborate a brand-new design method for electronically controlled bit-type units. First, based on the basic principle of the eigenmode theory in scattering, by analogy with the relationship between the total induced current and the eigencurrent and the relationship between the total scattered field and the eigenfield, it can be known that the phase of the eigenfield can represent the phase of the reflection coefficient, and the technical principle of designing an inverting digitalized unit based on the eigenmode theory parameters (mode expansion coefficient, eigenangle, polarization angle) is derived. Through this digitalized phase condition, the designer's expected physical characteristics of the automated design result are more clearly defined, and the electromagnetic characteristics of the final result with this as the optimization objective function are more in line with the design expectations. Combining this technical principle with the genetic optimization mechanism, first, the grid processing method in the preset area of the unit metal patch layer is used to complete the operations of retaining and removing the structure, and a certain number of structural individuals with different shapes are obtained as the initial population. Then, based on the MATLAB language, the algorithm is written to calculate the eigenmode parameters. According to the digitalized phase condition characterized by the eigenmode parameters, the automated design of digitalized units is directly completed through iteration, and at the same time, the purpose of simplifying the structure design process and shortening the optimization process is achieved.

[0061] The technical solution of the present invention describes an automated design method for a digital metasurface unit based on physical parameters of characteristic mode theory, comprising the following steps:

[0062] Step 1: Preset the basic structure of the unit;

[0063] A reflective digital metasurface unit structure loaded with a PIN diode, characterized in that it is a basic structural unit consisting of a first metal layer 1, a first dielectric layer 2 and a second metal layer 3. Figure 2 , the metal layer refers to the first metal layer 1 and the second metal layer 3 distributed from top to bottom; the dielectric layer is the first dielectric layer 2. According to the center frequency f of the expected working frequency band of the metasurface unit, the corresponding wavelength at this frequency is λ, and the unit structure period P can be preset to range as follows: Select the PIN diode model. In the present invention, the first metal layer 1 is a metal patch layer, the second metal layer 3 is a full-coverage metal layer under the unit structure period, and the first dielectric layer 2 is a dielectric substrate material.

[0064] Step 2: Preset the optimization area of ​​genetic iteration and mesh the metal patch layer;

[0065] The step 2 takes into account the PIN diode model selected in step 1, and presets the position of adding the diode. The present invention adopts the method of loading the diode by slit in the center of the first metal layer 1, and sets the position of introducing the diode to the center of the structure; then presets the desired design form of the basic structure of the unit. The present invention adopts a multi-resonant slit structure to ensure the working bandwidth of the unit. Therefore, a rectangular metal sheet containing a "艹"-shaped slit is selected as the initial structural form of the first metal layer 1. By gridding the metal patch layer, its initial structure is cut into the form of countless identical small grids of Nx rows and Ny columns, where the "艹"-shaped slit is the optimized area in the first metal layer 1 that needs to be randomly retained and deleted. The present invention adopts the form of dividing into small square grids, and the grid side length ml should be less than one twentieth of the wavelength λ at the center frequency f. The grid size is as small as possible, which will make the genetic optimization algorithm contain more gene bits, thereby improving the structural accuracy to ensure the individual diversity within the population. The first metal layer 1 is gridded in the present invention to facilitate the subsequent operation of step 3. The small grids within the optimization area preset in step 2 are randomly coded and assigned 0 or 1 to indicate whether the small grid is retained. A value of 0 indicates that the small grid is deleted, and a value of 1 indicates that the small grid is retained. The grids outside the preset optimization area of ​​the first metal layer 1 are always coded and assigned 1, that is, the grids are always retained, thereby achieving the purpose of completing structural optimization only within the predetermined optimization area.

[0066] Step 3: Set the number of individuals in the genetic algorithm population POP and the number of iterations GEN;

[0067] In Step 3, the number of individuals in the initial population is set to POP. By generating POP random binary arrays, an initial population consisting of POP unit structure grid individuals with different structures is obtained. The number of iterations is set to GEN. The larger the number of individuals POP in the initial population, the more individuals can be selected from the parents. The larger the number of iterations GEN, the closer the final optimization accuracy is to the target. However, too many population individuals or a large number of iterations will increase the computational load and extend the process of automated design. Appropriate numbers of individuals and iterations need to be set to complete the automated design faster.

[0068] Step 4: Set the objective function of the algorithm.

[0069] The research object for which the objective function set in Step 4 needs to be calculated is the unit structure form represented by each grid individual in the population obtained after Steps 1, 2, and 3. The characteristics of the objective function are related to the electromagnetic physical properties (characteristic mode phase difference) of the desired 1-bit unit structure. The objective function set in Step 4 is to complete the calculation of the objective function for the structure form represented by each grid individual in the population in Step 5. Therefore, the specific settings of the two objective functions for genetic optimization in Step 4 are as follows:

[0070] Objective function function 1: It is specified that the polarization mode of the incident wave is Y polarization or X polarization. It is required that the characteristic main mode polarization of the metasurface unit structure grid in the on-off state is the same as the polarization mode of the incident wave, and the target polarization angle polar should tend to 90° or -90°, that is, as shown in Equation (1):

[0071] function 1 = |polar| - 90 (1)

[0072] Objective function function 2: Since the digital phase difference of the characteristic mode is determined by the superposition of the phase difference of the main mode expansion coefficients, the difference in characteristic angles, and the difference in polarization angles of the two states of the electronic control unit under the same incident wave, the digital phase condition of the characteristic mode of the electronic control digital unit with a desired reflection field phase difference of Δgoal phase is:

[0073] function 2 = ΔMWC + ΔCA + Δpolar - Δgoal phase (2)

[0074] As described in Equation (2), where ΔMWC is the phase difference of the main mode pattern expansion coefficients in the on and off states, ΔCA is the difference between the main mode characteristic angles of the two, Δpolar is the difference between the main mode polarization angles of the two, and Δgoalphase is the desired phase difference of the design target. The minimum value of the optimization target optimizes to make function 2 tend to 0. With iterative optimization, the characteristic mode digital phase difference characteristic of the target digital unit will approach the desired phase difference of the design target.

[0075] Step 5: Solve the characteristic modes of the grid structure individuals representing the metal patch layer structure in the population and calculate the objective function;

[0076] The above Step 5 is based on the expected optimization objective function represented by the relevant characteristic mode physical parameters (main mode pattern expansion coefficient, characteristic angle, polarization angle) that can reflect the characteristic field phase in Step 4. Combining with the NSGA-ii genetic optimization algorithm, the script is written using the MATLAB language, and the model establishment, parameter solution, and data result recording work are completed by calling the electromagnetic simulation software CST Studio Suite 2020;

[0077] 5-1) Model establishment

[0078] Generate POP random binary arrays as described in Step 3. Each array is in the form of a discrete grid of an individual structure represented by random binary information. The POP grid individuals form the initial population. Sequentially read the binary information of the POP grid individuals and call the electromagnetic software CST Studio Suite 2020 to complete the construction of the structure model at the corresponding positions of the grid;

[0079] First, call the CST software and create a new project file. Considering that each grid individual in the initial population is a random binary row vector array, this row vector array is converted into an array matrix with Nx rows (total number of grid rows) and Ny columns (total number of grid columns). Then read this array matrix row by row and judge whether the grid in the optimization area is retained according to the binary information (the grids outside the optimization area are all retained, as described in Step 2). If retained, construct the grid structure model at the corresponding position until the construction of the overall upper metal grid structure corresponding to the entire array matrix is completed. Then add a new medium according to the parameters of the intermediate dielectric layer board and set the dielectric constant, and complete the construction of the intermediate dielectric layer structure according to the unit period P and the dielectric thickness h; finally, complete the construction of the fully covered metal floor structure on the lower surface of the dielectric layer according to the unit period size. Thus, the construction of the metasurface unit structure based on one grid individual in the population is completed;

[0080] 5-2) Parameter solution

[0081] For the metasurface units represented by each grid individual, the copy and translation operations are completed according to the unit period P. After arranging them into 4×4 finite periodic structures, the field monitors at the frequency f are set: far field farfield, magnetic field h-field, and then the engineering solver is set to a high-frequency integral equation solver (this engineering solver supports characteristic mode solving, and sets the number of characteristic mode modes ModeNum for solution analysis). Then, the setting of the Y-polarized plane wave excitation source incident vertically along the -Z axis to the unit structure is completed, and the electric field vector 'EVector' is set to (0,1,0). If the expected plane wave excitation source is an X-polarized wave incident vertically, just change the electric field vector 'EVector' to (1,0,0). In addition, considering that the unit structure is a reflective metasurface unit and the lower surface of the unit is a fully covered metal floor, the engineering model boundary condition Zmin is set to 'electric'. At this point, the solution setting of the project is completed, and the engineering solution of the characteristic mode is started.

[0082] 5-3) Data result recording

[0083] When the project completes the solution calculation, it is necessary to record the characteristic mode parameter data used for the objective function calculation, including the amplitude and phase of the mode expansion coefficient, the characteristic angle, and the characteristic field. This process is achieved by calling CST with MATLAB to read the project results and record the data. The far field, characteristic angle, and amplitude and phase data of the mode expansion coefficient of the characteristic mode are recorded one by one according to the characteristic mode sequence number. The data results are exported in ASCII format according to the project path, named in sequence according to the mode sequence number, and saved in each parameter folder.

[0084] At this point, the unit structure model is established, the characteristic mode parameters are solved, and the data result recording is completed, and the characteristic mode solution of the structural individuals in the population in step 5 is completed.

[0085] In order to realize the calculation of the objective function in step 5, it is necessary to determine the main mode of the characteristic mode and complete the extraction of relevant data of the objective function. The main mode of the characteristic mode has two characteristics: the mode expansion coefficient is relatively large; the polarization mode of the characteristic mode is the same as the polarization mode of the excitation source. Based on the above two characteristics, the specific process is as follows:

[0086] First, read the mode expansion coefficient MWC saved in the previous step, find the characteristic mode ordinal number with the largest amplitude of the mode expansion coefficient, and read and record the phase of the mode expansion coefficient of this mode ordinal number as MWC_phase. Then, read the X-component amplitude |Ex| and Y-component amplitude |Ey| in the far-field component data of this mode ordinal number. If it is excited by a Y-polarized wave, the mode polarization angle is denoted as: polar = tan(|Ey| / |Ex|), and for excitation by an X-polarized wave, the mode polarization angle is denoted as: polar = tan(|Ex| / |Ey|). Then, read the characteristic angle of this mode ordinal number and denote it as CA. For each grid individual in the population in step three, there are two structural forms under the on-off state of the diode. The only difference between the two is whether there is a grid structure at the diode position in the engineering structure model. Therefore, both can use the above steps to complete the solution and recording of the characteristic mode parameters in the two on-off states of the diode. Thus, the polarization angle polar for objective function 1, the phase difference ΔMWC of the mode expansion coefficient, the difference ΔCA of the characteristic angle, and the difference Δpolar of the polarization angle for objective function 2 can be calculated, so as to complete the calculation of the two objective functions. After recording the objective function results in the binary gene bit information of the grid individual, it is convenient for step six to select the best to generate the offspring population according to the objective function values in the subsequent steps.

[0087] Step six: Generate offspring population individuals through genetic iteration;

[0088] In step six, according to the objective function values recorded in the individual gene bit information, non-dominated sorting and crowding degree calculation are completed. The specific process is as follows:

[0089] 6-1) Non-dominated sorting

[0090] According to the objective function values, the individuals in the population are hierarchically prioritized through non-dominated sorting. The lower the non-dominated level number to which an individual belongs, the lower the priority and the higher the dominant position. The specific process of non-dominated sorting is as follows:

[0091] (1) The population size is N, and the objective function values of the individuals in the population are x(1), x(2) …… x(N);

[0092] (2) Let i = 1;

[0093] (3) For all j = 1, 2 …… N and j ≠ i, judge the dominance relationship between individual i and individual j by comparing the objective function x(i) of the i-th individual and the objective function x(j) of the j-th individual. If there is no x(j) that is better than x(i), then x(i) is denoted as a non-dominated individual;

[0094] (4) Let i = i + 1, and repeat to (2) until all non-dominated individuals are found.

[0095] The non-dominated individual set obtained through the above processes (1) to (4) serves as the first non-dominated layer of the population (with a priority of 0). Then, the individuals in the population other than the first non-dominated layer are repeatedly subjected to (1) to (4) to obtain the second non-dominated layer (with a priority of 1), and so on until the dominance relationship of the entire population is hierarchically sorted.

[0096] 6-2) Crowding Degree Calculation

[0097] Calculating the crowding degree is to estimate the crowding degree of other solutions around a specific individual solution in the population. The crowding degree is the average distance of adjacent individual solutions on both sides calculated according to each objective function. The specific calculation process of the crowding degree is as follows:

[0098] (1) Sort the population individuals in ascending order according to the magnitude of each objective function value;

[0099] (2) For each objective function, perform non-dominated sorting on the population. Take the individual solutions corresponding to the maximum value fmax and the minimum value fmin of the objective function as the boundary solutions, and set the crowding degree of the boundary individual solutions to an infinite distance. The crowding degree of other intermediate individual solutions is the normalized absolute difference of the function values of the adjacent individual solutions on both sides. Denote the m-th objective function value of the j-th individual solution as f (m) (j). In the case where this individual solution does not belong to the boundary solution, its crowding degree value is expressed as (f (m) (j + 1) - f (m) (j - 1)) / (f (m) max - f (m) min);

[0100] After non-dominated sorting and crowding degree calculation, each individual in the population has these two attributes: non-dominated layer (i.e., priority) and crowding degree value, and this is used as the basis for comparing the advantages and disadvantages of individuals in the population;

[0101] If the non-dominated layer to which the i-th individual x(i) belongs is superior to the non-dominated layer to which the j-th individual x(j) belongs, then the individual x(i) is superior to x(j); if x(i) and x(j) belong to the same non-dominated layer, then it is necessary to compare their crowding degrees. If the crowding degree f (m) (i) of the m-th objective function of x(i) is greater than the crowding degree f (m) (j) of the m-th objective function of x(j), then x(i) is superior to x(j). Therefore, according to the non-dominated sequence size and crowding degree of individuals, select individuals preferentially to form a new population for subsequent iteration, and repeat step five for the individuals in the new population until the objective function is optimal or until the maximum number of iterations is reached and stop.

[0102] Step Seven: Obtain the optimal individual, and the structural automatic design process ends.

[0103] The process of the genetic optimization algorithm in Step 7 ends. The on-off two states corresponding to the optimal individual are the digital unit that can be automatically designed and used as a reflective phase shifter. In Step 2 mentioned above, the position where the PIN diode is added is determined. Just introduce the equivalent circuit of the diode's conduction or truncation directly here. Thus, the automatic design of the encodable metasurface unit loaded with diodes is completed based on the eigenmode phase condition. The electromagnetic characteristics of the digital unit formed after the optimal structure designed automatically is in the form of the equivalent circuit of the loaded diode can meet the eigenmode digital phase at the central operating frequency f, meeting the design conditions.

[0104] The present invention first derives and proposes the digital phase condition based on the eigenmode theory, and summarizes that the eigenmode digital phase condition for designing the electrically controlled digital unit is jointly determined by the superposition of three parts: the phase difference of the main mode expansion coefficients, the difference in eigen angles, and the difference in polarization angles under two states of the electrically controlled unit under the same incident wave. As the technical principle for designing the metasurface unit loaded with active electrically controlled components, it no longer analyzes the characteristics of the reflection coefficient phase in the traditional "black box" manner. Compared with the existing design of electrically controlled 1-bit units, the method proposed by the present invention designs from the physical characteristics, and the developed optimization design process based on this can better reflect the physical characteristics of the structure.

[0105] Compared with the existing design methods, the technical method of the present invention combines the digital phase condition characterized by the physical concept parameters of eigenmodes as the objective function of the design, and uses the genetic iterative optimization algorithm mechanism as the technical means to directly complete the structure design based on the structure form, phase, frequency band, etc. conditions set by the designer according to the expected target function. During the process, there is no need for manual structure adjustment and optimization. The technical solution of the present invention not only makes the electromagnetic physical characteristics of the final optimization target more in line with the expected design target, but also greatly reduces the workload of the designer and shortens the design process.

[0106] Embodiment 1

[0107] Refer to Figure 1 , combined with Embodiment 1, the specific steps of the technology of the present invention are described in detail.

[0108] Step 1: Preset the basic structure of the unit

[0109] Refer to Figure 2 , a reflective digital metasurface unit structure loaded with a PIN diode, characterized in that it is a basic structural unit composed of three parts: a first metal layer 1, a first dielectric layer 2, and a second metal layer 3. The first metal layer 1 and the second metal layer 3 are distributed from top to bottom; the dielectric layer is the first dielectric layer 2. The first metal layer 1 and the second metal layer 3 are respectively located on the upper surface and the lower surface of the first dielectric layer 2.

[0110] The unit structure period P is 12 mm, which is 0.4 times the wavelength λ. In this embodiment, the first dielectric layer 2 is made of a dielectric material with a thickness h of 1.58 mm and a relative dielectric constant ε r The first metal layer 1 is a Taconic TLX substrate material with a loss tangent tanδ of 2.55 and a loss tangent tanδ of 0.0019. Figure 3 The side length L of the square metal patch is 7.5 mm, the center of the patch is slotted and loaded with a PIN diode of Skyworks SMP1340. When the diode is turned on, its equivalent circuit is a series circuit of a 1Ω resistor and a 450pH inductor; when it is cut off, its equivalent circuit is a series circuit of a 10Ω resistor, a 450pH inductor and a 100fF capacitor. The second metal layer 3 is a metal plate that fully covers the unit cycle.

[0111] Step 2: Preset the “艹”-shaped optimization region for genetic iteration

[0112] Reference Figure 3 , the square metal patch of the first metal layer 1 is divided into small grid areas with Nx being 30 rows and Ny being 30 columns, and the side length ml of each small grid is 0.25 mm. Considering the unit polarization mode and diode model, the "艹"-shaped optimization area inside the first metal layer 1 is composed of gap 1, gap 2 and gap 3.

[0113] The gap 1 is a preset gap along the X-axis in the 15th row in the center of the patch grid area, with a width ml of 0.25mm and a length L of 7.5mm, wherein the PIN tube is located at the small grid position of the 15th row and the 15th column. When the PIN tube is turned on and off, the electronic control unit can be ideally equivalent to the presence or absence of the metal sheet, so the code assignment of 0 or 1 at the small grid of the 15th row and the 15th column can represent on and off. In order to ensure that the final unit structure has a relatively wide working bandwidth, the slits are symmetrically opened along the Y axis to increase the resonant structure. The gap 2 on the left side of the Y axis is characterized by a length L of 7.5 mm and a width d2 of 1.5 mm, corresponding to the 9th to 14th columns of the grid area, and a distance d1 of 2 mm from the left edge of the metal sheet; the gap 3 on the right side of the Y axis is characterized by a length L of 7.5 mm and a width d2 of 1.5 mm, corresponding to the 16th to 21st columns of the grid area, and a distance d1 of 2 mm from the right edge of the metal sheet, and together with the gap 1 (the 15th row of the grid area) it forms a "艹"-shaped gap. This "艹"-shaped gap area (except the grid point where the PIN tube is located) is defined as the optimization area for genetic iteration.

[0114] Step 3: Set the number of individuals in the genetic algorithm population POP and the number of iterations GEN;

[0115] The small grids in the "艹"-shaped optimization area in step 2 are randomly coded with 0 or 1 to indicate whether the small grid is retained (0 means the small grid is deleted, 1 means the small grid is retained), while the grids outside the optimization area of ​​the first metal layer 1 are always coded with 1 (grid retention), thereby achieving the purpose of structural optimization only in the "艹"-shaped area, generating many unit grids with different structures as individuals in the initial population of the algorithm.

[0116] In Example 1, the number of individuals in the initial population POP is set to 40, and an initial population consisting of 40 unit structure grid individuals with different structures is obtained; the number of iterations GEN is set to 20. The more the number of initial population POP, the more individuals that can be selected from the parents, and the larger the number of iterations GEN, the closer the final optimization accuracy is to the target. However, if there are too many individuals in the population or the number of iterations is large, the amount of calculation will increase and the process of automated design will be prolonged. Therefore, it is necessary to set an appropriate number of individuals and iterations in order to complete the automated design faster.

[0117] Step 4: Set the objective function of the algorithm;

[0118] Objective function function 1: The polarization mode of the incident wave is specified as Y polarization, and the characteristic main mode polarization corresponding to the metasurface unit structure grid in the on-off state is required to be the same polarization mode as the incident wave. The target polarization angle polar should tend to 90° or -90°, as shown in formula (1):

[0119] function 1=|polar|-90 (1)

[0120] Objective function 2: According to the characteristic mode digitization phase difference, which is determined by the phase difference of the main mode mode expansion coefficient, the difference of the characteristic angle and the difference of the polarization angle under the same incident wave, the characteristic mode digitization phase condition of the electrically controlled 1-bit unit with a reflection field phase difference Δgoal phase of 180° is:

[0121] function 2=ΔMWC+ΔCA+Δpolar-Δgoal phase (2)

[0122] As described in formula (2), ΔMWC is the phase difference of the main mode expansion coefficient in the on and off states, ΔCA is the difference between the main mode characteristic angles, and Δpolar is the difference between the main mode polarization angles. The minimum optimization algorithm is adopted in the technical solution of the present invention, and the optimization target function 2 tends to 0. In Example 1, it is stipulated that the characteristic mode phase difference of 1bit unit characteristics is required to be 180°, that is, the Δgoal phase of the expected phase difference representing the design target is 180°. According to the designer's expected phase requirements, it can be modified according to needs.

[0123] Step 5: Solve the characteristic modes of the grid structure individuals representing the metal patch layer structure in the population, and calculate the objective function;

[0124] In this embodiment, Step 5 is based on the expected optimization objective function represented by the relevant characteristic mode physical parameters (main mode expansion coefficient, characteristic angle, polarization angle) that can reflect the phase of the characteristic field in Step 4. Combining with the NSGA-ii genetic optimization algorithm, the script is written using the MATLAB language. The model establishment, parameter solution, and data result recording work are completed by calling the electromagnetic simulation software CST Studio Suite 2020. The specific process is as follows:

[0125] 5-1) Model establishment

[0126] Generate POP random binary arrays according to Step 3. Each array is a discrete grid form of an individual structure represented by random binary information. The POP grid individuals form the initial population. Sequentially read the binary information of the POP grid individuals, and call the electromagnetic software CST Studio Suite 2020 to complete the construction of the structural model at the corresponding positions of the grids.

[0127] First, call the CST software and create a new project file. Considering that each grid individual in the initial population is a random binary row vector array, this row vector array is converted into an array matrix with 30 rows (total number of grid rows) Nx and 30 columns (total number of grid columns) Ny. Then read this array matrix row by row, and judge whether the grids in the optimization area are retained according to the binary information (the grids outside the optimization area are all retained, as described in Step 2). If retained, construct the grid structure model at the corresponding positions until the construction of the overall upper metal grid structure corresponding to the entire array matrix is completed. Then add a new medium according to the parameters of the intermediate dielectric layer board and set the dielectric constant, and complete the construction of the intermediate dielectric layer structure according to the unit period P and the dielectric thickness h; finally, complete the construction of the fully covered metal floor structure on the lower surface of the dielectric layer according to the unit period size. Thus, the construction of the metasurface unit structure based on one grid individual in the population is completed.

[0128] 5-2) Parameter solution

[0129] For the metasurface units represented by each grid individual, the copy and translation operation is completed according to the unit period P, and after arranging them into 4×4 finite periodic structures, the field monitors at the frequency f are set: far field (f=10GHz), magnetic field h-field (f=10GHz), and then the engineering solver is set to a high-frequency integral equation solver (this engineering solver supports characteristic mode solving, and the number of characteristic mode modes ModeNum for solution analysis is set to 20), and then the setting of the Y-polarized plane wave excitation source incident vertically along the -Z axis is completed. In addition, considering that the unit structure designed by the present invention is a reflective metasurface unit, and the lower surface of the unit is a fully covered metal floor, the engineering model boundary condition Zmin is set to 'electric'. At this point, the solution setting for the engineering is completed, and the engineering solution of the characteristic mode begins.

[0130] 5-3) Data result recording

[0131] When the project completes the solution calculation, it is necessary to record the characteristic mode parameter data used for the objective function calculation, including the amplitude and phase of the mode expansion coefficient, characteristic angle, and characteristic field. This process is achieved by calling CST in MATLAB to read the project results and record the data. The far field, characteristic angle, and amplitude and phase data of the mode expansion coefficient of the characteristic mode are recorded one by one according to the characteristic mode sequence number, and the data results are exported in ASCII format according to the project path, and named in sequence according to the mode sequence number and saved in each parameter folder.

[0132] At this point, the unit structure model is established, the characteristic mode parameters are solved, and the data result recording is completed, and the characteristic mode solution of the structural individuals in the population in step 5 is completed.

[0133] To calculate the objective function in Step 5, it is necessary to determine the dominant mode of the eigenmode and complete the extraction of relevant data for the objective function. The dominant mode of the eigenmode has two characteristics: the modal expansion coefficient is relatively large; the polarization mode of the eigenmode is the same as that of the excitation source. Based on the above two characteristics as the judgment basis, the specific process is as follows: First, read the modal expansion coefficient MWC saved in the previous step, find the eigenmode ordinal number with the largest amplitude of the modal expansion coefficient, and read and record the phase of the modal expansion coefficient of this mode ordinal number as MWC_phase; then read the X-component amplitude |Ex| and Y-component amplitude |Ey| in the far-field component data of this mode ordinal number, and the modal polarization angle is recorded as: polar = tan(|Ey| / |Ex|); then read the characteristic angle of this mode ordinal number and record it as CA. For each grid individual in the population in Step 3, there are two structural forms under the on-off state of the diode, and the difference between the two is only whether there is a grid structure at the diode position in the engineering structure model. Therefore, the two can respectively use the above steps to complete the solution and recording of the eigenmode parameters in the two on-off states of the diode. Therefore, the polarization angle polar for the objective function 1, the phase difference ΔMWC of the modal expansion coefficient, the difference ΔCA of the characteristic angle, and the difference Δpolar of the polarization angle for the objective function 2 can be calculated, so as to complete the calculation of the two objective functions. After recording the objective function results in the binary gene position information of the grid individual, it is convenient for the subsequent Step 6 to select the best according to the objective function value to generate the offspring population.

[0134] Step 6: Generate offspring population individuals through genetic iteration;

[0135] Perform non-dominated sorting and crowding degree calculation on the individuals in the population. According to the non-dominated sequence size and crowding degree of the individuals, select the best to generate a new population for subsequent iteration, and repeat Step 5 for the individuals in the new population until the objective function is optimal or the maximum number of iterations is reached and stop. The above Step 6 completes non-dominated sorting and crowding degree calculation according to the objective function values recorded in the individual gene position information. The specific process is as follows:

[0136] 6-1) Non-dominated sorting

[0137] Perform priority stratification on the individuals in the population through non-dominated sorting according to the objective function values. The lower the non-dominated level number to which an individual belongs, the lower the priority and the higher the dominant position. The specific process of non-dominated sorting is as follows:

[0138] (1) The population size is N, and the objective function values of the individuals in the population are x(1), x(2)... x(N);

[0139] (2) Let i = 1;

[0140] (3) For all j = 1, 2,..., N and j ≠ i, the dominance relationship between individual i and individual j is determined by comparing the objective function x(i) of the i-th individual and the objective function x(j) of the j-th individual. If there is no x(j) that is superior to x(i), then x(i) is recorded as a non-dominated individual;

[0141] (4) Let i = i + 1, and repeat step (2) until all non-dominated individuals are found.

[0142] The set of non-dominated individuals obtained through the above processes (1) to (4) serves as the first non-dominated layer (priority 0) of the population. Then, the other individuals in the population except for the first non-dominated layer are repeatedly processed through (1) to (4) to obtain the second non-dominated layer (priority 1), and so on until the dominance relationship of the entire population is hierarchically sorted.

[0143] 6-2) Crowding degree calculation

[0144] Calculating the crowding degree is to estimate the crowding degree of other solutions around a specific individual solution in the population. The crowding degree is the average distance of adjacent individual solutions on both sides calculated according to each objective function. The specific process of calculating the crowding degree is as follows:

[0145] (1) Sort the population individuals in ascending order according to the magnitude of each objective function value;

[0146] (2) For each objective function, perform non-dominated sorting on the population. Take the individual solutions corresponding to the maximum value fmax and the minimum value fmin of the objective function as the boundary solutions, and set the crowding degree of the boundary individual solutions to an infinite distance. The crowding degree of other intermediate individual solutions is the normalized absolute difference of the function values of the adjacent individual solutions on both sides. Denote the m-th objective function value of the j-th individual solution as f (m) (j). In the case where this individual solution does not belong to the boundary solution, its crowding degree value is expressed as (f (m) (j + 1) - f (m) (j - 1)) / (f (m) max - f (m) min).

[0147] After non-dominated sorting and crowding degree calculation, each individual in the population has these two attributes: non-dominated layer (i.e., priority) and crowding degree value, which are used as the basis for comparing the quality of individuals in the population. If the non-dominated layer to which the i-th individual x(i) belongs is superior to the non-dominated layer to which the j-th individual x(j) belongs, then individual x(i) is superior to x(j); if x(i) and x(j) belong to the same non-dominated layer, then it is necessary to compare their crowding degrees. If the crowding degree f (m) (i) of the m-th objective function of x(i) is greater than the crowding degree f (m)(j), then x(i) is superior to x(j). Therefore, according to the non-dominated sorting sequence size and crowding degree of individuals, individuals are preferentially selected to form a new population for subsequent iterations. Step five is repeated for the individuals in the new population until the objective function is optimal or the maximum number of iterations is reached and then stopped.

[0148] Step seven: Obtain the optimal individual, and the structural automatic design process ends.

[0149] Refer to Figure 4 、 Figure 5 illustrates the metal patch structure area of the electronic control unit corresponding to the grid of the optimal individual. Step two clarifies the addition position of the PIN diode, so the equivalent circuit of the conduction or truncation of the diode can be directly introduced here. Therefore, the on and off states corresponding to the optimal individual obtained after the iteration termination can be obtained, which is the 1-bit digital unit that can be used as a reflective phase shifter automatically designed. Refer to Figure 4 is the grid structure of the first metal layer 1 in the conduction state. Inside a square metal sheet with a side length L of 7.5 mm, there are two gaps along the Y-axis with a length L1 of 6.5 mm and a width S of 0.5 mm (the distance from the gap edge is L2 of 2 mm), and a gap along the X-axis in the center of the patch with a square grid representing the position of the PIN tube blocking in the middle. Figure 5 is the grid structure of the first metal layer 1 in the cut-off state, with two gaps along the Y-axis with a length L1 of 6.5 mm and a width S of 0.5 mm (the distance from the gap edge is L2 of 2 mm), and a gap along the X-axis in the center with a length L2 of 2 mm and a width ml of 0.25 mm. Thus, the automatic design of the encodable metasurface unit loaded with diodes is completed based on the eigenmode phase condition. The electromagnetic characteristics of the digital unit formed after the optimal structure automatically designed is in the form of the equivalent circuit of the loaded diode can meet the eigenmode digital phase at the center operating frequency f of 10 GHz, meeting the design conditions.

[0150] The parameters used in this embodiment are shown in Table 1.

[0151] Table 1 Parameter list in Embodiment 1

[0152] P L d1 d2 12 mm 7.5 mm 2 mm 1.5 mm ml L1 L2 S 0.25 mm 6.5 mm 2 mm 0.5 mm

[0153] In Embodiment 1 of the present invention, the preset frequency point f is 10 GHz, and the characteristic mode phase difference is 180°. However, the embodiments based on the technical method of the present invention do not necessarily adopt the setting conditions of a frequency point of 10 GHz and a characteristic mode phase difference of 180°. The unit period P in Step 1, the shape of the optimization region in Step 2, the objective function in Step 4, etc. can be modified according to the desired frequency band and phase effect of the designer to meet the design requirements. Embodiment 1 of the present invention does not necessarily adopt a unit period P of 12 mm. Embodiment 1 of the present invention does not necessarily adopt the set population individual number POP and iteration number GEN, which can be modified according to the actual iteration process to meet the design requirements. In Embodiment 1 of the present invention, it is not necessary to adopt the Taconic TLX substrate material with the same parameters and the same model PIN diode. After considering costs and processing, etc., replacements can be made according to the actual application requirements such as the board material, components, and application frequency band to meet the design requirements.

Claims

1. An automated design method for digital metasurface units based on characteristic mode theory, characterized in that Including the following steps; Step 1: Preset the basic structure of the digital metasurface unit; Step 2: Adopt the genetic algorithm, preset the optimization region of genetic iteration according to the basic structure obtained in Step 1, and perform grid processing on the metal patch layer; Step 3: Set the population individual number POP and the iteration number GEN of the genetic algorithm for the grid-processed metal patch layer; Step 4: Correspondingly set the objective function of the algorithm for the unit structure represented by each grid individual in the population obtained in Step 3; Step 5: Solve the characteristic modes of the grid structures representing the metal patch layer structures in the population, and calculate the objective function; Step 6: Generate offspring population individuals through genetic iteration of the objective function obtained in Step 5; Step 7: Obtain the optimal individual, and the automated structure design process ends; In order to calculate the objective function in Step 5, it is necessary to judge the main mode of the characteristic mode and complete the extraction of relevant data for the objective function. The main mode of the characteristic mode has two characteristics: the mode expansion coefficient is relatively large; the polarization mode of the characteristic mode is the same as the polarization mode of the excitation source. Based on the above two characteristics as the judgment basis, the specific process is as follows: First, read the mode expansion coefficient MWC saved in the previous step, find the characteristic mode ordinal number with the largest amplitude of the mode expansion coefficient, and read and record the phase of the mode expansion coefficient of this mode ordinal number as MWC_phase; then read the X-component amplitude |Ex| and Y-component amplitude |Ey| in the far-field component data of this mode ordinal number, and the mode polarization angle is recorded as: polar = tan(|Ey| / |Ex|); then read the characteristic angle of this mode ordinal number and record it as CA. For each grid individual in the population in Step 3, there are two structural forms under the on-off state of the diode. The difference between the two is only the presence or absence of the grid structure at the diode position in the engineering structure model. Therefore, the two can respectively use the above steps to complete the solution and recording of the characteristic mode parameters in the two on-off states of the diode. Therefore, the polarization angle polar for the objective function 1, the phase difference ΔMWC of the mode expansion coefficient, the difference ΔCA of the characteristic angle, and the difference Δpolar of the polarization angle for the objective function 2 can be calculated, so as to complete the calculation of the two objective functions. After recording the objective function results in the binary gene bit information of the grid individual, it is convenient for the subsequent Step 6 to select the best according to the objective function value to generate the offspring population.

2. The automated design method of the digital metasurface unit based on the characteristic mode theory according to claim 1, characterized in that The basic structure of the first step is a reflective digital metasurface unit structure loaded with PIN diodes, which is a basic structural unit composed of three parts: a first metal layer (1), a first dielectric layer (2), and a second metal layer (3). According to the center frequency f of the desired operating frequency band of the metasurface unit structure, with the corresponding wavelength being λ, the preset range of the unit structure period P is and the PIN diode model is selected.

3. The automated design method of the digital metasurface unit based on the characteristic mode theory according to claim 1, characterized in that In the step 2, the position for adding the diode is preset in consideration of the PIN diode model selected in the step 1, and the diode introduction position is set to the center of the structure by means of a slot loading method in the center of the first metal layer (1); then the desired design form of the basic structure of the unit is preset, and a multi-resonant slot structure is adopted to ensure the working bandwidth of the unit, and a rectangular metal sheet containing a "艹"-shaped slot is adopted as the initial structure form of the first metal layer (1); and the initial structure is cut into the form of countless identical small grids of Nx rows and Ny columns by means of a gridding method of the metal patch layer, wherein the "艹"-shaped slot is The first metal layer (1) is an optimized area where random mesh retention and deletion are required. The area is divided into small square meshes. The mesh side length ml should be less than one twentieth of the wavelength λ at the center frequency point f. The small meshes in the optimized area preset in step 2 are randomly coded and assigned 0 or 1 to indicate whether the small mesh is retained. The assigned value 0 indicates that the small mesh is deleted, and the assigned value 1 indicates that the small mesh is retained. The meshes outside the preset optimized area of ​​the first metal layer (1) are always coded and assigned 1, that is, the meshes are always retained, thereby achieving the purpose of completing the structural optimization only in the predetermined optimized area.

4. The automated design method of the digital metasurface unit based on the characteristic mode theory according to claim 1, wherein In the step three, the number of individuals in the initial population is set to POP, and an initial population consisting of POP unit structure grid individuals with different structures is obtained by generating POP random binary arrays; the number of iterations is set to GEN, and the more the number of POP in the initial population means the more individuals that can be selected from the parents, and the larger the number of iterations GEN means the closer the final optimization accuracy is to the target.

5. The automated design method of the digital metasurface unit based on the characteristic mode theory according to claim 1, characterized in that, The research object for calculating the objective function in step 4 is the unit structure form corresponding to each grid individual in the population obtained after steps 1, 2, and 3, and the characteristics of the objective function are related to the electromagnetic physical properties of the final desired digital unit structure. The objective function set in step 4 is to complete the objective function calculation under the structural form represented by each grid individual in the population in step 5. The specific settings of the two objective functions of genetic optimization in step 4 are as follows: Objective function function 1: The polarization mode of the incident wave is specified as Y polarization or X polarization, and the characteristic main mode polarization corresponding to the metasurface unit structure grid in the on-off state is required to be the same as the incident wave. The target polarization angle polar should tend to 90° or -90°, as shown in formula (1): function1=|polar|-90 (1) Objective function function 2: According to the characteristic mode digitization phase difference is determined by the superposition of the main mode mode expansion coefficient phase difference, characteristic angle difference and polarization angle difference of the two states of the electronically controlled unit under the same incident wave, the characteristic mode digitization phase condition of the electronically controlled digital unit with the reflection field phase difference of Δgoal phase is: function2=ΔMWC+ΔCA+Δpolar-Δgoal phase (2) As described in Equation (2), where ΔMWC is the phase difference of the main mode pattern expansion coefficients in the on and off states, ΔCA is the difference in the main mode characteristic angles between the two, Δpolar is the difference in the main mode polarization angles between the two, and Δgoal phase is the desired phase difference of the design goal. The optimization process of minimizing the optimization target makes function2 tend to 0. With iterative optimization, the characteristic mode digital phase difference characteristic of the target digital unit will approach the desired phase difference of the design goal.

6. The automated design method of the digital metasurface unit based on the characteristic mode theory according to claim 1, characterized in that Step 5 is based on the desired optimization objective function represented by the relevant characteristic mode physical parameters that can reflect the phase of the characteristic field in Step 4. Combining with the NSGA-ii genetic optimization algorithm, the script is written using MATLAB language, and the model establishment, parameter solution, and data result recording work are completed by calling the electromagnetic simulation software CST Studio Suite 2020. 5-1) Model establishment Generate POP random binary arrays according to Step 3. Each array is in the form of a discrete grid of an individual structure represented by random binary information. The POP grid individuals form the initial population. Sequentially read the binary information of the POP grid individuals and call the electromagnetic software CST Studio Suite 2020 to complete the construction of the structural model at the corresponding positions of the grid. First, call the CST software and create a new project file. Considering that each grid individual in the initial population is a random binary row vector array, this row vector array is converted into an array matrix with Nx rows and Ny columns. Then read the array matrix row by row and determine whether the grid in the optimization area is retained according to the binary information. If it is retained, construct the grid structure model at the corresponding position until the construction of the overall upper metal grid structure corresponding to the entire array matrix is completed. Then add a new medium according to the parameters of the intermediate dielectric layer plate and set the dielectric constant, and complete the construction of the intermediate dielectric layer structure according to the unit period P and the dielectric thickness h. Finally, complete the construction of the fully covered metal floor structure on the lower surface of the dielectric layer according to the unit period size. Thus, the construction of the metasurface unit structure based on one grid individual in the population is completed. 5-2) Parameter solution For the metasurface unit represented by each grid individual, complete the replication and translation operations according to the unit period P, arrange it into a 4×4 finite periodic structure, and set the field monitors at frequency f: far field farfield and magnetic field h-field. Then use the high-frequency integral equation solver as the engineering solver, and complete the setting of the Y-polarized plane wave excitation source incident vertically on the unit structure along the -Z axis. Set the electric field vector ’EVector’ to (0,1,0). If the desired plane wave excitation source is an X-polarized wave incident vertically, just change the electric field vector ’EVector’ to (1,0,0). In addition, considering that the unit structure is a reflective metasurface unit and the lower surface of the unit is a fully covered metal floor, set the engineering model boundary condition Zmin to ’electric’. Thus, the solution setting of the project is completed, and the engineering is started to solve the characteristic modes. 5-3) Data result recording After the project completes the solution calculation, it is necessary to record the eigenmode parameter data used for the objective function calculation, including the amplitude and phase of the mode expansion coefficients, eigenangles, and eigenfields. This process is achieved by MATLAB calling CST to read the project results and record the data. The farfield, eigenangles, and the amplitude and phase data of the mode expansion coefficients of the eigenmodes are recorded successively according to the eigenmode order. The data results are exported in ASCII format according to the project path and named successively according to the mode order, and stored in each parameter folder. So far, the establishment of the unit structure model, the solution of the eigenmode parameters, and the recording of the data results have been completed, and the eigenmode solution of the structural individuals in the population in Step 5 has been completed.

7. The automated design method of the digital metasurface unit based on the characteristic mode theory according to claim 1, wherein In Step 6, according to the objective function values recorded in the individual gene position information, non-dominated sorting and crowding degree calculation are completed. The specific process is as follows: 6-1) Non-dominated sorting Through non-dominated sorting based on the objective function values, the individuals in the population are stratified by priority. The lower the non-dominated level number to which an individual belongs, the lower the priority and the higher the dominant status. The specific process of non-dominated sorting is as follows: (1) The population size is N, and the objective function values of the individuals in the population are x(1), x(2)... x(N); (2) Let i = 1; (3) For all j = 1, 2... N and j ≠ i, the dominance relationship between individual i and individual j is judged by comparing the objective function x(i) of the i-th individual and the objective function x(j) of the j-th individual. If there is no x(j) that is better than x(i), then x(i) is recorded as a non-dominated individual; (4) Let i = i + 1, and repeat to (2) until all non-dominated individuals are found; The set of non-dominated individuals obtained through the above processes (1) to (4) is used as the first non-dominated layer of the population (priority 0). Then, the other individuals in the population except the first non-dominated layer are repeated for (1) to (4) to obtain the second non-dominated layer (priority 1), and so on until the dominance relationship of the entire population is sorted hierarchically; 6-2) Crowding degree calculation Calculating the crowding degree is to estimate the crowding degree of other solutions around a specific individual solution in the population. The crowding degree is the average distance between adjacent individual solutions calculated according to each objective function. The specific process of crowding degree calculation is as follows: (1) Sort the population individuals in ascending order according to the magnitude of each objective function value; (2) For each objective function, perform non - dominated sorting on the population. Take the individual solutions corresponding to the maximum value fmax and the minimum value fmin of the objective function as the boundary solutions. Set the crowding distance of the boundary individual solutions to infinite distance, and set the crowding distance of other intermediate individual solutions to the normalized absolute difference of the function values of the two adjacent individual solutions on both sides. Denote the m - th objective function value of the j - th individual solution as f (m) (j). In the case where this individual solution does not belong to the boundary solution, its crowding distance value is expressed as (f (m) (j + 1)-f (m) (j - 1)) / (f (m) max - f (m) min); After non-dominated sorting and crowding degree calculation, each individual in the population has these two attributes: non-dominated layer (i.e., priority) and crowding degree value, which are used as the basis for comparing the quality of individuals in the population. If the non-dominated layer to which the \(i\)-th individual \(x(i)\) belongs is superior to the non-dominated layer to which the \(j\)-th individual \(x(j)\) belongs, then the individual \(x(i)\) is superior to \(x(j)\); if \(x(i)\) and \(x(j)\) belong to the same non-dominated layer, then it is necessary to compare their crowding degrees. If the crowding degree \(f\) (m) (i) of the \(m\)-th objective function of \(x(i)\) is greater than the crowding degree \(f\) (m) (j) of the \(m\)-th objective function of \(x(j)\), then \(x(i)\) is superior to \(x(j)\). Therefore, according to the size of the non-dominated sequence and the crowding degree of the individuals, individuals are preferentially selected to form a new population for subsequent iteration. The steps in Step 5 are repeated for the individuals in the new population until the objective function is optimal or until the maximum number of iterations is reached and then stopped.

8. The automated design method of the digital metasurface unit based on the characteristic mode theory according to claim 1, characterized in that In Step 7, the optimal individual is obtained, and the genetic optimization algorithm process ends. The on-off two states corresponding to the obtained optimal individual are the digital units that can be automatically designed and used as reflective phase shifters. In Step 2, the addition positions of the PIN diodes are specified. Just introduce the equivalent circuit of the diode's conduction or truncation directly here. Thus, the automated design of the encodable metasurface unit loaded with diodes is completed based on the eigenmode phase condition. The electromagnetic characteristics of the digital unit formed after the optimal structure automatically designed is in the form of the equivalent circuit of the loaded diode can satisfy the eigenmode digital phase at the central operating frequency f and meet the design conditions.

Images