• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    A Numerical Convex Lens for the State-Discretized Modeling and Simulation of Megawatt Power Electronics Systems as Generalized Hybrid Systems

    2021-04-22 11:35:04BochenShiZhengmingZhaoYichengZhuZhujunYuJiaheJu
    Engineering 2021年12期

    Bochen Shi, Zhengming Zhao*, Yicheng Zhu, Zhujun Yu, Jiahe Ju

    Department of Electrical Engineering, Tsinghua University, Beijing 100084, China

    Keywords:Generalized hybrid systems Megawatt power electronics Modeling and simulation Numerical convex lens

    ABSTRACT Modeling and simulation have emerged as an indispensable approach to create numerical experiment platforms and study engineering systems.However,the increasingly complicated systems that engineers face today dramatically challenge state-of-the-art modeling and simulation approaches. Such complicated systems,which are composed of not only continuous states but also discrete events,and which contain complex dynamics across multiple timescales, are defined as generalized hybrid systems (GHSs) in this paper. As a representative GHS, megawatt power electronics (MPE) systems have been largely integrated into the modern power grid, but MPE simulation remains a bottleneck due to its unacceptable time cost and poor convergence.To address this challenge,this paper proposes the numerical convex lens approach to achieve state-discretized modeling and simulation of GHSs. This approach transforms conventional time-discretized passive simulations designed for pure-continuous systems into state-discretized selective simulations designed for GHSs.When this approach was applied to a largescale MPE-based renewable energy system, a 1000-fold increase in simulation speed was achieved, in comparison with existing software.Furthermore,the proposed approach uniquely enables the switching transient simulation of a largescale megawatt system with high accuracy, compared with experimental results,and with no convergence concerns.The numerical convex lens approach leads to the highly efficient simulation of intricate GHSs across multiple timescales, and thus significantly extends engineers’capability to study systems with numerical experiments.

    * Corresponding author.

    E-mail address: zhaozm@mail.tsinghua.edu.cn (Z. Zhao).

    https://doi.org/10.1016/j.eng.2021.07.011

    2095-8099/? 2021 THE AUTHORS. Published by Elsevier LTD on behalf of Chinese Academy of Engineering and Higher Education Press Limited Company.This is an open access article under the CC BY-NC-ND license(http://creativecommons.org/licenses/by-nc-nd/4.0/).

    1. Introduction

    With the rapid development of computational technologies in recent decades,numerical modeling and simulation have emerged as an indispensable approach to recreate,explore,decode,and predict complex dynamics of interest across a wide range of scientific disciplines and engineering applications [1,2]. However, with the rapid development of science, technology, and modern society,the dynamic behavior of the physical systems of interest to humanity is becoming increasingly complicated, such that these systems can no longer be simply regarded as pure-continuous systems [3,4]. The emergence of such complicated systems dramatically challenges state-of-the-art modeling and simulation approaches.

    Fig.1 provides some examples.In ecosystems(Fig.1(a)),events such as crustal movement, volcanic eruption, and glacial melting fundamentally change the evolution mode of the system, particularly when taking into account climate change due to recent human activities [5,6]. Another example is the spread and control of coronavirus disease 2019 (COVID-19) as shown in Fig. 1(b), in which control measures such as travel bans and city shutdowns changed from time to time according to the spread of infection,and in turn governed the spread pattern of the disease [7,8]. Engineering examples include chemical processing plants, in which continuous states such as temperature interact with discrete events such as phase transition[9].In robot control,foot impacting and knee locking lead to mode transitions [10]. Clearly, the above systems cannot simply be modeled as pure-continuous systems.They show a hybrid nature involving the coexistence of both continuous states and discrete events.

    Fig. 1. Typical examples of generalized hybrid systems (GHSs) and their multi-timescale nature. (a) Discrete events in an ecosystem fundamentally change the evolution mode(e.g.,biodiversity).The timescale of these events is relatively very small.(b)Control measures for COVID-19 as discrete events change the infection modes of the virus as continuous states.Peak mobility may occur during the short period of the discrete event(e.g.,after the announcement and before its execution).‘‘Week/div”refers to the unit of the x-axis is one week per grid. (c) In power converters, switching events of the semiconductor switches control the transfer of energy flow as desired. The nanosecond-level switching transient causes an abrupt change of electromagnetic energy and, sometimes, system failure, especially in megawatt power electronics (MPE)systems.

    Another engineering example that we focus on in this paper is that of power electronics systems (Fig. 1(c)). Power electronics have been widely used in nearly all electric energy conversion fields [11]. In particular, given the increasing concerns about climate change and the energy crisis, megawatt power electronics(MPE) systems act as indispensable interfaces to connect renewable sources with the utility grid and end users[12].In power electronics converters, the originally continuous voltage and current are‘‘discretized”through the introduction of power semiconductor switches [13]. Hence, the system demonstrates the coexistence of both continuous states and discrete events, and can thus be defined as a hybrid system.

    The term ‘‘hybrid system” originated from a special workshop that was held in 1986 [14], and soon became an emerging topic[4,15,16].The early defined hybrid system mainly involved continuous states controlled by discrete signals [17]. Thus far, the concept has been enormously expanded to include the mixing of and interaction between continuous states and discrete events[15]. However, the abovementioned definitions and discussions still cannot depict the full complexity of the hybrid nature of engineering systems today.On many occasions,discrete events cannot occur instantaneously.It is known that time is the measurement of existence; therefore, there must be a time interval during the‘‘event,”although it is often too fast to be perceived.From a physical perspective, a discrete event still has a continuous transient with a small timescale; this results in a prominent multitimescale feature in generalized hybrid systems (GHSs) (Fig. 1).The transients of discrete events can have great significance,especially when the physical energy during the short transient of an event is large [18]. A typical example is the switching transient in MPE (Fig. 1(c)). Because the power level is high, the ‘‘sudden change”and‘‘imbalance”in large electromagnetic energy is highly likely to cause system faults [19–21].

    These switching transients cannot be described within the conventional concept of a hybrid system and the ideal concept of a discrete event without any physical process. Therefore, in this paper,we introduce and define the concept of a GHS in order to describe a more general and complete context of the hybrid nature of physical systems. We define a GHS as follows: A GHS is a multitimescale dynamic system composed of both time-driven continuous dynamics and event-driven discrete dynamics, where largetimescale continuous dynamics embody discrete processes, and discrete processes embody small-timescale continuous transients.

    In GHSs, continuous states and discrete events as two distinct types of variables co-exist and interact with each other, collaboratively determining the system behavior. Meanwhile, two types of continuous states of different timescales simultaneously exist within the system, resulting in a prominent multi-timescale feature that must be considered in high-power applications. Such GHS are not only commonly found in the physical world, but also regarded as a more general form of the dynamic system.Subsequently,a challenge of such a complicated dynamic behavior—which is of equal significance to many engineers in various applications—is the effective and efficient modeling and simulation of GHSs.

    The remainder of this paper is organized as follows: The challenges of modeling and simulating GHSs using conventional approaches are summarized in Section 2. A so-called ‘‘numerical convex lens” designed based on the features of a GHS in order to achieve state-discretized modeling and simulation is then proposed in Section 3 to address these challenges. The proposed approach has general value for different GHSs in various engineering applications. More specifically, we examine how to model and simulate an MPE system as a representative GHS. Applications of the proposed approach in MPE systems are demonstrated in Section 4. Finally, conclusions are drawn in Section 5.

    2. Challenges of modeling and simulating GHSs

    2.1. Modeling and simulation challenges for general GHSs

    In current modeling and simulation approaches [22], time is taken as the independent variable to build the system model and determine the corresponding solver; this is called the timediscretized approach. This approach is a natural idea: Since time is the measurement of existence, to study a dynamic system is to explore the pattern of the system dynamics based on time. However, time flows continuously and uniformly without qualitative change. In contrast, the variables of state in GHSs can not only accumulate, but also change qualitatively. In other words, when a certain condition is fulfilled, the system will switch to another mode and continue its accumulation under a totally different pattern. In this situation, numerical experiments based on timediscretization encounter a significant challenge: The discretized timepoints may not be the exact time when an event occurs.Therefore,such experiments are incapable of evaluating the event’s impact on the system dynamics accurately. As a consequence, an extremely small timestep must be adopted in order to accurately locate the discrete event and solve the corresponding smalltimescale transient during the event,which results in dramatically low efficiency, or even computational failure.

    2.2. Challenges of modeling and simulation for MPE systems

    The subsection above provided a general discussion of the challenges encountered when modeling and simulating GHSs. We now discuss MPE systems in detail in order to demonstrate more concrete examples.The modeling of MPE systems can be basically categorized into three types according to different levels of semiconductor models:average models [23,24], ideal models [25–27], and physical models[28–30]. The average model eliminates (averages) the switching actions within one or several switching cycles and therefore transforms the GHS into a pure-continuous system. Thus, the analyses and designs of the system—especially those of the closed-loop controls—can be well-supported with available mathematical tools such as the transfer function and the bode plot [31]. Within the context of GHSs,the average model only focuses on the large-timescale continuous dynamics. The impact of discrete events is ignored, let alone switching transients within the events.When considering the switching events, the ideal switch model is widely used in many systemlevel simulation tools [26]. This model considers the ‘‘largetimescale continuous–discrete”process,but still ignores the switching transients of switching events.

    Existing models of switching transients are complicated equivalent circuits based on semiconductor physics, as shown in Fig. 2[28]. Such physical models abandon the assumption of discrete events and describe the semiconductor switches as smalltimescale continuous systems instead. The GHS is then transformed again into a pure-continuous system, but this time into a small-timescale system that considers the switching transients.Theoretically, physical models accurately describe smalltimescale transients. Practically, these physical models have also been implemented in device-level commercial tools for power electronics [29,32] and have been used in studies with a small number of switches. Unfortunately, it is hardly possible to use them in MPE applications, in which typically hundreds or thousands of semiconductor devices exist and interact with each other[33,34].An example is the simulation of a 24-switch 50 kVA solidstate transformer with physical switch models tested in Ref. [35];although it is just a medium-scale circuit, the simulation already takes about 9 h and frequently encounters divergence failures during numerical studies.It is scarcely possible to simulate larger converters in MPE simulations with physical models within an acceptable time period while ensuring a convergent result.

    To sum up, it is still challenging to model the integral‘‘continuous–discrete–continuous” process in MPE. Existing approaches are still based on time as the variable to develop the models; therefore, only a certain aspect (i.e., continuous/discrete)and certain timescales can be focused on.

    The solving of MPE systems as GHSs is equally challenging due to the enormous number of switching events.Existing approaches such as state-space methods [26,27] and nodal analysis methods [36],which are equipped with various fixed- or variable-step integration methods[37],still fall under the category of time-discretized solvers.Such methods are incapable of automatically matching discretized points with the occurrence of discrete events and,furthermore,with the small-timescale transients within the events, which are fundamental features introduced by GHSs.Iterations are frequently necessary during the simulation in order to locate the discrete event, and small timesteps must be adopted to solve the transients and overcome convergence issues, as shown in Fig. 3. For MPE systems with hundreds or thousands of semiconductor switches, the switching events of different switches are typically asynchronous. Even for a relatively low switching frequency, such as 1 kHz, 1000 switches in the simulated system can cause events to occur 1 million times per second. This results in low simulation speed and poor convergence.

    3. A numerical convex lens for state-discretized modeling and simulation

    Fig. 2. Existing modeling approaches for MPE systems include the average model, ideal model, and physical model. Each model focuses on one aspect of the ‘‘continuous–discrete–continuous” process of the GHS,and the models vary from large to small timescales. For the average model,D is the duty cycle, and L,C, and R are the inductance,capacitance,and resistance in a buck converter,respectively.For the physical model,Cgs,Cgd,and Cdsj are the metal-oxide-semiconductor capacitances,Imult,Icss,and Ibss are the equivalent current sources in insulate gate bipolar transistor(IGBT),Ccer and Ceb are the equivalent capacitances in IGBT,Rb is the equivalent resistance connected to the anode, and Imos is the metal-oxide-semiconductor current.

    To address the above challenges of modeling and simulation for GHSs, we propose the concept of state-discretized modeling and simulation.The general idea involves role-swapping the variable of time and the variable of state. In pure-continuous systems, time can be a sufficient variable to perform modeling and simulation.However,in GHSs,time is only a passive variable.Due to the introduction of discrete events and the corresponding small-timescale transients within the events, system states become a more dominating variable, which not only embodies the system dynamics,but also triggers the occurrence of events. By swapping the roles of time and state, a state-discretized solver can potentially merge the discretization of the simulation with the quantitative change of the GHS(i.e.,the occurrence of the events),and is therefore capable of modeling the integral process of multiple timescales while maximizing the simulation efficiency.

    Fig. 3. Existing time-discretized solvers for MPE systems. Discretized points in a time-discretized solver can be mismatched with the occurrence of the discrete events.Iterative calculations are necessary to locate the events,and small timesteps must be adopted to simulate switching transients without convergence problems.

    In this section, in order to implement the concept of statediscretized modeling and simulation, we receive inspiration from optics and propose a novel method called the ‘‘numerical convex lens.” Although they are considered to be two separate fields,optics and computational science share the same target: to perform observation,either with light or with computations.In optical experiments,if a certain part of an object cannot be clearly seen,a convex lens can be used to magnify it in detail.The degree to which the convex lens magnifies the observed object depends not only on the lens itself,but also on the incident rays from the object,such as their light frequency and intensity. This analogy inspired us to compare the system variables in numerical experiments(including continuous states and discrete events) to light rays and endow them with different ‘‘frequency” and ‘‘intensity,” so that they will have different refractive indexes through the numerical convex lens (i.e., the solver) and, therefore, different levels of magnification. In this way, we can pick out variables that have undergone qualitative changes (which means that their magnified images are clear enough to be observed) and let the clear imaging trigger the numerical experiment. This leads to a ‘‘selective” numerical observation, where we swap the roles of state and time. To be specific, with the help of the convex lens, the discretization of states(i.e.,the imaging process)triggers the solving of the system,and the discretization of time is determined accordingly. In this way, we can efficiently match the discretized solving points with the qualitative changes of system states without unnecessary calculations,which enables accurate and fast numerical experiments.

    3.1. A fractal model of the GHS with the numerical convex lens

    To perform numerical observation, the first task is numerical modeling. With the assistance of the numerical convex lens, we can establish a ‘‘fractal model” for GHSs (Fig. 4). In the fractal model, system dynamics are decoupled into different groups according to their differences in timescale, which are defined as different dynamic planes. The dynamics in the fast planes are the transient processes of discrete events in slow planes.In each plane,system dynamics can be modeled with the hybrid automaton model [25,38]. When observing the fractal model using the numerical convex lens, the system variables will be imaged through the lens. This process is called ‘‘state discretization.” The clear image triggers the numerical solving. In particular, when a discrete event occurs, the numerical convex lens will magnify it into a new dynamic plane and switch the numerical experiment into this plane, where the small-timescale behavior of the system is described as a new ‘‘small hybrid system.” When this discrete event ends, the numerical experiment will be switched back to the previous plane. In this way, the dynamics of different timescales can be decoupled,so that the solving speed and convergence of the numerical experiment can be simultaneously improved with high accuracy.The above idea is very similar to the idea of fractals in geometry:Research on fractal geometry has revealed that when observing the geometric features of a system on different spacescales, it is usually found that parts exhibit an amazing similarity to the whole; likewise, in numerical experiments, with the help of the numerical convex lens,we can observe the system dynamics in different timescales on different dynamic planes, and the dynamics on different planes can be similarly modeled as hybrid systems. Therefore, we call this numerical model of multitimescale GHSs the ‘‘fractal model.”

    In the fractal model,system dynamics are decoupled into different groups according to their differences in timescale and defined as different dynamic planes. For example, as shown in Fig. 4, the system dynamics are divided into three groups and defined as three different dynamic planes: ①the slow dynamic plane containing second-to millisecond-scale dynamics;②the fast dynamic plane containing microsecond-to nanosecond-scale dynamics;and③ the ultra-fast dynamic plane containing picosecond- to femtosecond-scale dynamics. To observe the details in a discrete event in the slower plane, the numerical convex lens can be used to magnify the discrete event into a faster plane,which represents a small hybrid system. In other words, the dynamics in the fast planes are the switching transient of discrete events in the slow planes. If necessary, we can also magnify the switching transients in this small hybrid system into a dynamic plane with a ‘‘smaller timescale.”By repeating the above process,we can use the numerical convex lens to observe the dynamic behavior of the generalized dynamic system in different timescales on different dynamic planes.In this way,dynamics with different timescales can be distinguished, thus avoiding the need to consider unnecessary dynamic behaviors.

    In each plane, the hybrid system can be modeled as a hybrid automaton [38], as illustrated in Fig. 5. The hybrid automaton model comprises the following elements:

    (1) Variables

    ①Time variable t.

    ②State variables: A set of variables inside the system that can describe the dynamic behavior of the system completely and without redundancy.

    Fig. 4. Establishing the fractal model of GHSs based on the numerical convex lens.GHSs exhibit prominent multi-timescale features. Based on the numerical convex lens, we can decouple the system dynamics into different groups and define them as different dynamic planes. Discrete events in each plane can be magnified and modeled as a new dynamic plane by the numerical convex lens.

    Fig. 5. The hybrid automaton model for a single dynamic plane. The model comprises four elements: variable, mode, guard, and event. A dynamic plane contains a set of all possible dynamic modes and guard conditions.

    · Continuous state variables x: State variables that vary continuously with time.

    · Discrete state variables d: State variables that vary discontinuously with time.

    ③Input variables u:Variables that are outside the system but can influence the dynamic behavior inside the system. Input variables can be continuous or discrete.

    (2)Mode:The operating pattern of systems dynamics.A combination D of discrete variables dictates a mode. In different modes,the dynamic patterns are different, which means that the corresponding mathematical models are different. For example, in Fig. 5, for the kth mode mk, we can describe the system mode as the following ordinary differential equation (ODE):

    where c is called the characteristic variable and can be expressed by state variables, and vthis the corresponding threshold value. When Eq. (2) holds true, system states will remain in the present mode.When the characteristic variable crosses the threshold,it will cause a change of the inequality’s Boolean value, and thus may trigger a mode transition of the hybrid system.

    (4) Event: The mode transition of the hybrid system. An event happens when the guard condition of the present mode changes its Boolean value.

    A dynamic plane contains a set of all possible dynamic modes and guard conditions.When a discrete event occurs and its transition details are of concern, the numerical convex lens can be utilized to magnify it and transfer the numerical experiment to the faster plane. When the transition process ends, an ending event will be triggered,and the numerical experiment can be transferred back to the slower plane.

    3.2.State-discretized solving of a GHS with the numerical convex lens

    Fig. 6. Observing the fractal model of GHSs through the numerical convex lens.When solving the fractal model, the color(frequency) of the continuous states and the light intensity of the discrete events reflect the degree of attention that should be paid in the numerical experiment.They exhibit different refractive indices when they pass through the numerical convex lens,such that variables with a clear image,which are the variables that undergo qualitative changes, can be picked out to spontaneously trigger the numerical experiment. yblue and yred are the image heights of the blue and red light,respectively,and fblue and fred are the focal lengths of the blue and red light, respectively.

    This section introduces how to solve the fractal model of GHSs.According to the principles of the numerical convex lens (Fig. 6),the system variables will separate into different light rays; hence,clearly imaged variables that have undergone qualitative changes will be picked out and will spontaneously trigger numerical experiments. In general, to solve a GHS composed of continuous states and discrete events,different imaging principles are used for these two categories of elements with different properties. Below are detailed descriptions on how to solve GHSs with the numerical convex lens.

    Let u0be the object distance, y0be the object height, and y be the image height in Fig. 6. Let us assume that, in order to be captured by human eyes,the height of the image should be more than the resolution threshold ylim, that is,

    Then the numerical criteria for clear image formation can be rewritten as

    The remaining question becomes: What is the definition of the refractive index n for system variables? As was mentioned in the introduction of this paper,a GHS comprises two types of system variables: continuous states and discrete events. Since these two elements have different properties, their imaging principles in the numerical convex lens system are different.In the numerical system,we compare them to rays of light and endow continuous states with different colors (frequencies) and discrete events with different intensities. According to the refraction principle of the numerical convex lens,lights with different colors or intensities will have different refractive indices and different focal lengths, and their images will naturally be separated in space and have different sizes.

    For continuous states,taking a single-variable case as an example,we compare it to a beam of monochromatic light with variable color,whose color at a certain moment corresponds to a frequency that is called the ‘‘numerical angular frequency.” Here, we define the numerical angular frequency as the measurement error in numerical experiments. At present, there have been many algorithms designed for pure-continuous systems that can provide the local error estimate at each timestep, such as the widely used ode45 algorithm[37]and others[39,40].Similar to the Cauchy dispersion formula in optics, the refractive index of the numerical convex lens for continuous states Δn1has the following nonlinear relationship with the numerical angular frequency ω:

    where cbaseis the base value for the characteristic variable c. Eq.(10) shows that the closer the characteristic variable of a discrete event is to the threshold, the greater the weight of the discrete event will be, and the higher the priority it should be given in numerical experiments. Real lens materials can have nonlinearity and exhibit different refractive indexes for light with different intensities,in what is called the Kerr effect in optics.Inspired by this phenomenon, we endow the numerical convex lens with a similar effect and define the relationship between the refractive index Δn2and the light intensity as follows:

    where K is the Kerr coefficient of the lens.The larger K is,the stronger the Kerr effect of the lens will be.Thus,K reflects another property of the lens,which represents the error tolerance of the solving of discrete events.

    Taking the image formation principle of both continuous states and discrete events into account, we can obtain the complete expression of the refractive index, as follows:

    For continuous states,we consider them to be variable-frequency and constant-intensity, and let KI <<1; for discrete events, we consider them to be variable-intensity and constant-frequency,and let ω2/b2<<1. Based on the numerical criteria for clear image formation given by Eq.(8)and the complete expression of the refractive index given by Eq.(12),the fractal model can be solved according to the formed image of the system variables in the generalized hybrid system. Since the numerical convex lens can have both the dispersion effect and the Kerr effect, it has the ability to efficiently solve continuous states and discrete events at the same time.

    3.3. Modeling and solving MPE systems

    To apply the above-described general methods to the simulation of MPE systems,this subsection discusses detailed interpretations of all the concepts in the numerical convex lens (including the fractal model and the solving algorithm) based on a power electronics context. To establish the fractal model shown in Fig. 4 for MPE systems,we divide the system dynamics into two dynamic planes according to their timescales.The first-level dynamic plane describes the large-timescale system-level dynamics, where the switching-on and switching-off of the power switches are considered as discrete events.When the switch is in its steady state(onstate or off-state, modeled as small resistance or open-circuit,respectively),the system can be described with the state equation:

    where x is the state vector,and each element of x is an independent state variable, namely the voltage of a capacitor and the current of an inductor; and u is the input vector, and each element of u is an input of the system, such as a power source. Moreover, it is shown in Eq.(13)that the state equation fswis dependent on sw,a switching function vector where each element is the switching state(on or off) of a switch.

    For mode transition, namely, the switching event, the model must define the guard for each switching event,as Eq.(2)suggests.Here, we list the main guard conditions for switching events.

    For a diode switch, the guard condition for its turn-on is when its voltage where vDis greater than zero:

    The guard condition for an active switch(e.g.,an insulated gate bipolar transistor, or IGBT) is determined by the controller of the system.In power electronics,pulse-width modulation(PWM)control is usually used to determine the switching signals[41],where switching events are triggered when a reference waveform vref(which depends on system states) crosses a carrier waveform vcar(a regular signal). Thus, the guard condition of an active switch is

    where vrefand vcarare the values of the reference and the carrier,respectively.

    Thus far, the first-level dynamic plane in the fractal model has been established. To describe the transient of a discrete event in the first-level dynamic plane—namely, the switching transient of the semiconductor device—it is necessary to build a hybrid model of the second-level dynamic plane. Here, in order to model the switching behavior, we use the piecewise analytical transient(PAT)modeling method as described in Ref.[20].In the PAT model,the switching transient is divided into different stages according to the physical mechanisms inside the semiconductor device(Fig. 7).Each stage corresponds to one mode (one circle in Fig. 5) in the second-level dynamic plane.

    Now that we have established the fractal model for MPE systems,with level-one and level-two dynamic planes,the next target is to solve the fractal model with the numerical convex lens. Following the principles of the numerical convex lens, our tasks include: ①to determine the numerical angular frequency for the continuous state variables; ②to determine the numerical light intensity for the discrete event variables; and ③to solve the system following the imaging principles. Below are descriptions of how these concepts are defined in MPE systems.

    First, as defined above, the numerical angular frequency is the numerical measurement error and is dependent on the integration algorithm. For non-stiff circuits of MPE systems, we can use the flexible adaptive discrete state (FA-DS) integration algorithm to calculate the numerical solution of the continuous states [39].According to the algorithm, the numerical angular frequency is defined as follows:

    Fig.7. PAT model of an IGBT-diode switching pair.The complicated physical mechanisms within semiconductors are decoupled,as shown in the left circuit.Only dominant mechanisms are considered in each stage of the switching transient,as shown in the right waveform.The relationship between the mechanisms and the stages are marked with colors. BJT refers to the bipolar junction transistor. Ls is the stray inductor and vLs is the voltage across it.CDC is the direct current-bus capacitor and VDC is the voltage across it.VGon,VGoff,RGon,and RGoff are the gate-drive voltages and resistances.RGint is the inner resistance of the gate node.Cgc,Cge,and Cce are the junction capacitances of the IGBT and vcg,vge,and vce are the voltages across them.DS is the diode,vD is the voltage across it,irr is the reverse recovery current of it,and iD is the current through it.Lload and Rload are the inductance and resistance of the load,and IL is the current through the load.ic is the current through the IGBT.ig is the gate-dive current of the IGBT.Icmax is the maximum current of the IGBT during the turn-on transient. t0 ~t6 are the time instants of each stage.

    where x is one of the state variables (the inductor current of the capacitor voltage), xbaseis the base value of x, p is the order of the integration algorithm,xp+1denotes the(p+1)th-order derivatives of x,and Δt denotes how long the numerical integration lasts since the last calculation timepoint.

    Based on the above definitions,ω(x,t)represents the estimated numerical error of the integration algorithm for state variable x at timepoint t,thus offering what the numerical convex lens requires for the imaging of continuous states: namely, the numerical measurement error(which denotes how much this state variable matters and whether we need to update the system states at this timepoint).

    For the discrete events, the definition of the numerical light intensity of each discrete event has been given in Eq. (10). We define the base value for the discrete event solving cbaseas the maximum value of all the threshold values defined in the fractal model—that is, Eqs. (14–16)—and the absolute tolerance as a setting of the algorithm.

    We can now solve the system following the image formation of the numerical convex lens, according to Eqs. (8–12) and Eq.(17). There are three parameters to control the simulation accuracy—that is, b, K, and β—which control the error tolerances of the continuous states, the discrete events, and the overall simulation, respectively. Once a clear imaging triggers a new calculation, we can use the discrete state algorithm to calculate the continuous states, or use the secant method [39] to locate the time of occurrence of discrete events. Furthermore, discrete events will be magnified into a second-level dynamic plane to simulate the switching transients.

    4. Applications and validations in MPE systems

    This section presents applications of the state-discretized approach in MPE systems. A 2 MW renewable energy system is taken as an example to test and verify the proposed modeling and simulation approach. The structure in Fig. 8 illustrates one of the operational modes of the system, with renewable energy and storage. The main component is a four-port electric energy router(EER)[34]that connects the utility grid,the renewables,the energy storage, and various types of loads.

    The topology of the EER is shown in Fig. 9. In the studied system,one EER contains 576 switching devices;in theory,this means that the system has 2576possible operating modes in total. In addition, each mode transition is accompanied by a nanosecondlevel electromagnetic transient. This makes the system unsolvable with any time-discretized approaches and tools currently available, unless the mathematical model is greatly simplified, which would result in low accuracy.Even with simplifications from ignoring the switching transients of the discrete events, the simulation still takes much too long (up to hours or days). Consequently,due to huge challenges in numerical modeling and solving, the design and analysis of such a system have mainly relied on experience or simplification in the past.

    To better elaborate the concept of the fractal model in the studied case,we show the mathematical form of the fractal model here,taking one high-voltage direct current (HVDC) submodule (Fig. 9)in the EER as an example.For the upper dynamic plane,if we define mode U1 as (SHD1= SHD4= 1, SHD2= SHD3= 0), and mode U2 as(SHD1= SHD4= 1, SHD2= SHD3= 0), we can derive the upper-level state equations within the fractal model as follows:

    Fig. 8. Structure of the studied MPE system. IT refers to information technology, PV refers to photovoltaic, DC refers to direct current and AC refers to alternating current.EER: electric energy router.

    Fig. 9. Topology of one EER. It contains a 10 kV high-voltage alternating current (HVAC) port, a 10 kV high-voltage direct current (HVDC) port, a ±375 V low-voltage direct current(LVDC)port and a 380 V low-voltage alternating current(LVAC)port.A 20 kHz high-frequency alternating current(HF-AC)link is used as the AC bus inside the EER.One EER contains 87 submodules (SMs) and 72 high-frequency transformers (HFTs). SHD1, SHD2, SHD3, and SHD4 are the switching signals of the SiC MOSFETs in a HVDC SM.LHD and CHD are the inductance and capacitance in a HVDC SM.

    Fig. 10. Numerical experiments on the 2 MW system. (a) Diagram of the prototype and internal structure of the EER, and a comparison between the principles of timediscretization and state-discretization methods in a numerical experiment of power electronic hybrid systems.Comparisons of(b)simulated results and(c)simulated speed between the numerical convex lens and commercial simulation software based on time-discretization. (d) Comparison of switching-current simulated results between the numerical convex lens and commercial simulation software that uses an ideal switch model. ns/div: nanosecond per division; μs/div: microsecond per division;ms: millisecond per division.

    where Lsis the stray inductance as in Fig. 7, and vSHD2,vSHD3, iSHD1,and iSHD4are the equivalent voltage and current sources of the MOSFETs.

    The expressions of vSHD2,vSHD3,vSHD1,and vSHD4are given by the PAT model[20].It can be found that the model for the upper-level dynamics—that is,Eqs.(18)and(19)—and the model for the lowerlevel transients—that is, Eqs. (20) and (21)—share very similar forms. This finding proves the concept of the fractal model; that is,when we zoom into the smaller timescale(Fig.4,the transients of the discrete events), the model shows similarity to the largertimescale model. These different timescales are all described by hybrid models,with multiple modes and states and algebraic equations in each mode.

    Based on the fractal model above, a numerical platform can be built for this EER system.A 0.2 s load-change dynamic of one EER is tested.Compared with the time-discretization approach,the statediscretization approach enabled by the numerical convex lens can dramatically improve the calculation efficiency (Fig. 10(a)). A professional simulation software in the field of power electronics commercial is chosen as the benchmark.All the tests are performed on the same personal computer, with a 4.20 GHz Intel Core i7-7700K CPU and 32 GB of memory. When simulating a 200 ms dynamic process (where switching transients are ignored because no existing tools can solve the system otherwise), the numerical convex lens achieves more than 1000 times acceleration at the same accuracy,shortening the simulation time from nearly 4 h to around 10 s(Figs.10(b)and(c)).This result is achieved on an ordinary personal computer and does not involve any acceleration technology by multi-core computation in parallel. Furthermore, the numerical convex lens takes advantage of the fractal model to solve the switching transient of discrete events—that is, the nanosecondlevel electromagnetic transient—thus attaining high consistency with experimental results (the simulation time with the proposed approach is then 608 s for a 0.2 s simulation when considering the switching transients of all the semiconductor switches).To date,no other commercial simulation software has such capability; other software can only approximately regard the switching events as instantaneous (Fig. 10(d)). The accuracy and high speed of the numerical convex lens enables a groundbreaking revolution in the design and development of the MPE system:All designers will be able to verify its design efficiently and accurately, solve multitimescale dynamics,and perform virtual control tests at low hardware cost (i.e., an ordinary personal computer) and low time cost(seconds or minutes). With the numerical convex lens, designers can conduct the previously unimaginable and time-consuming simulation of system dynamics, and can even achieve automatic design and iteration by introducing artificial intelligent algorithms.

    Another scenario of the EER is simulated to demonstrate the capability of the proposed numerical convex lens for addressing dynamic interactions across multiple timescales.In a real EER project, a high-frequency oscillation of 2 MHz is observed across the AC link[42].Due to the topology shown in Fig.10,this 2 MHz oscillation has a large influence on the stable operation of the system.It generates significant electromagnetic interference (EMI) and causes malfunctioning of the control system. Such highfrequency behavior is highly dependent on the small-timescale switching transients, while also having a large influence on the overall performance of the system. Therefore, this scenario provides typical evidence of the tight connection between multiple timescales in MPE systems [43]. With the numerical platform, we study the high-frequency oscillation with both a conventional gate driver(Fig.11(a))and an active gate driver(details of the design of the active driver can be found in Ref. [43]) (Fig. 11(b)). The simulated results show that with the active gate driver, which actively controls the voltage slew rate of the switching transient, it is largely possible to mitigate the high-frequency oscillations and therefore the corresponding EMI.

    5. Conclusions

    Although time-discretized simulation approaches have been well-developed for pure-continuous systems, the current bottleneck in the development of engineering systems lies in how to address emerging largescale complex systems, defined as GHSs,which involve both continuous states and discrete events across multiple time scales.To resolve this bottleneck,this paper demonstrated the numerical convex lens method to implement statediscretized modeling and simulation of GHSs. In the proposed approach,the calculated point is triggered by the imaging of a virtual convex lens, where the imaging process represents the discretization of the system states. Therefore, the discretized points can be naturally matched with the discrete events within GHSs.

    The proposed approach can be universally used for all GHSs among various engineering fields. In this work, we applied and verified the proposed approach in MPE systems,which are increasingly being integrated in modern power grids during the transition toward a more sustainable energy system. A 1000-fold speedup was achieved using the proposed approach. Furthermore, this approach has the unique capability to simulate nanosecond-level switching transients in such largescale systems without convergence problems. This approach will advance the computing capability and extend the limitations of numerical prototypes in engineering systems, such that simulations of complicated GHSs are no longer a bottleneck.

    Fig.11. Simulated results of high-frequency oscillation with(a)a conventional gate driver and (b) an active gate driver.

    Acknowledgments

    This work was supported by the Major Program of National Natural Science Foundation of China (51490683).

    Compliance with ethics guidelines

    Bochen Shi,Zhengming Zhao,Yicheng Zhu,Zhujun Yu,and Jiahe Ju declare that they have no conflict of interest or financial conflicts to disclose.

    999久久久精品免费观看国产| 国产亚洲欧美精品永久| 男女下面进入的视频免费午夜 | 午夜激情av网站| 9热在线视频观看99| 亚洲avbb在线观看| 无人区码免费观看不卡| 欧美在线一区亚洲| 亚洲国产欧美一区二区综合| 久久天躁狠狠躁夜夜2o2o| 久久久国产精品麻豆| 国产精品av久久久久免费| 日日干狠狠操夜夜爽| 丝袜美腿诱惑在线| 久久精品人人爽人人爽视色| 男人舔女人下体高潮全视频| 大型黄色视频在线免费观看| 亚洲人成网站在线播放欧美日韩| 黄色成人免费大全| 亚洲欧美日韩另类电影网站| 一边摸一边抽搐一进一出视频| 久热爱精品视频在线9| 亚洲国产精品合色在线| 999久久久精品免费观看国产| 亚洲成人精品中文字幕电影 | 99热只有精品国产| 黄色视频不卡| 18禁美女被吸乳视频| cao死你这个sao货| 亚洲色图 男人天堂 中文字幕| 午夜精品在线福利| 午夜精品久久久久久毛片777| 午夜精品在线福利| 在线观看免费午夜福利视频| 女警被强在线播放| 亚洲一区高清亚洲精品| 女人被狂操c到高潮| 中文字幕色久视频| 极品教师在线免费播放| 长腿黑丝高跟| 丝袜美腿诱惑在线| 国产色视频综合| 久久久精品国产亚洲av高清涩受| 国产区一区二久久| 色尼玛亚洲综合影院| 最新美女视频免费是黄的| 亚洲avbb在线观看| 国产精品 欧美亚洲| 最近最新中文字幕大全电影3 | 在线观看一区二区三区| 两个人免费观看高清视频| 99精品欧美一区二区三区四区| 90打野战视频偷拍视频| 久久精品国产亚洲av高清一级| 精品免费久久久久久久清纯| 1024香蕉在线观看| 99国产精品免费福利视频| 啦啦啦 在线观看视频| 巨乳人妻的诱惑在线观看| 19禁男女啪啪无遮挡网站| 国产91精品成人一区二区三区| 露出奶头的视频| 国产极品粉嫩免费观看在线| 热re99久久国产66热| 男人舔女人的私密视频| 国内毛片毛片毛片毛片毛片| 亚洲一区高清亚洲精品| 欧美乱码精品一区二区三区| 亚洲第一av免费看| 午夜福利在线观看吧| 亚洲男人的天堂狠狠| 曰老女人黄片| 国产av一区二区精品久久| 夜夜爽天天搞| 午夜成年电影在线免费观看| 黄色成人免费大全| 欧美午夜高清在线| 99在线视频只有这里精品首页| 色尼玛亚洲综合影院| 日韩大尺度精品在线看网址 | 色综合婷婷激情| 欧美日韩精品网址| 黄网站色视频无遮挡免费观看| 人人妻,人人澡人人爽秒播| 久久国产乱子伦精品免费另类| 精品一区二区三卡| 无遮挡黄片免费观看| 纯流量卡能插随身wifi吗| 久久欧美精品欧美久久欧美| 国产99白浆流出| 亚洲精品在线观看二区| 色在线成人网| 涩涩av久久男人的天堂| 精品电影一区二区在线| 日韩欧美国产一区二区入口| 欧美一区二区精品小视频在线| 精品电影一区二区在线| 999精品在线视频| 麻豆久久精品国产亚洲av | 十八禁网站免费在线| a级片在线免费高清观看视频| 亚洲欧美日韩无卡精品| 丝袜在线中文字幕| 国产乱人伦免费视频| 美女福利国产在线| 亚洲在线自拍视频| 老司机福利观看| 97超级碰碰碰精品色视频在线观看| 91国产中文字幕| 最近最新中文字幕大全免费视频| 国产三级黄色录像| 亚洲久久久国产精品| 日韩精品免费视频一区二区三区| 免费看a级黄色片| 纯流量卡能插随身wifi吗| 热re99久久精品国产66热6| 18禁观看日本| 亚洲国产毛片av蜜桃av| 精品一区二区三卡| 久久精品国产清高在天天线| 午夜a级毛片| 国产精品国产av在线观看| 男女下面进入的视频免费午夜 | 人成视频在线观看免费观看| 亚洲av熟女| 久久精品国产亚洲av高清一级| 亚洲国产欧美一区二区综合| 久久九九热精品免费| 如日韩欧美国产精品一区二区三区| 在线观看免费日韩欧美大片| tocl精华| 69精品国产乱码久久久| 中文字幕人妻熟女乱码| 成年人黄色毛片网站| 国产成人一区二区三区免费视频网站| e午夜精品久久久久久久| 真人一进一出gif抽搐免费| 91av网站免费观看| 日日干狠狠操夜夜爽| 麻豆国产av国片精品| 极品教师在线免费播放| 久久性视频一级片| 可以免费在线观看a视频的电影网站| 男男h啪啪无遮挡| 亚洲av五月六月丁香网| 欧美亚洲日本最大视频资源| 久久狼人影院| 女生性感内裤真人,穿戴方法视频| 国产三级黄色录像| 亚洲成av片中文字幕在线观看| 亚洲av日韩精品久久久久久密| 亚洲欧美一区二区三区黑人| 1024视频免费在线观看| 国产精品成人在线| 最近最新免费中文字幕在线| 免费看十八禁软件| 国产精品秋霞免费鲁丝片| 精品午夜福利视频在线观看一区| 久久天躁狠狠躁夜夜2o2o| 一进一出抽搐动态| 欧洲精品卡2卡3卡4卡5卡区| 久久久久精品国产欧美久久久| 午夜免费鲁丝| 丝袜美足系列| 国产精品乱码一区二三区的特点 | 正在播放国产对白刺激| 咕卡用的链子| 久久久久久人人人人人| 99久久人妻综合| 91成人精品电影| av视频免费观看在线观看| a级片在线免费高清观看视频| 亚洲成av片中文字幕在线观看| 天天添夜夜摸| 1024视频免费在线观看| tocl精华| 久久久精品国产亚洲av高清涩受| 首页视频小说图片口味搜索| 免费久久久久久久精品成人欧美视频| 亚洲情色 制服丝袜| 91在线观看av| 欧美不卡视频在线免费观看 | 久久久久国产精品人妻aⅴ院| 777久久人妻少妇嫩草av网站| 水蜜桃什么品种好| 色哟哟哟哟哟哟| 免费少妇av软件| 少妇被粗大的猛进出69影院| 欧美成人午夜精品| 侵犯人妻中文字幕一二三四区| 国产高清视频在线播放一区| 无遮挡黄片免费观看| 久久精品国产清高在天天线| 美女高潮到喷水免费观看| 精品国产乱码久久久久久男人| 久久久久久久精品吃奶| 亚洲精华国产精华精| 久久99一区二区三区| 亚洲黑人精品在线| 啪啪无遮挡十八禁网站| 一边摸一边抽搐一进一小说| 亚洲人成电影免费在线| 亚洲人成伊人成综合网2020| 日本wwww免费看| 免费久久久久久久精品成人欧美视频| 久久久久精品国产欧美久久久| 91九色精品人成在线观看| 亚洲中文字幕日韩| 自线自在国产av| 久热这里只有精品99| 国产黄a三级三级三级人| 亚洲七黄色美女视频| 午夜精品久久久久久毛片777| 法律面前人人平等表现在哪些方面| 久久久久久久久中文| 亚洲色图av天堂| 麻豆久久精品国产亚洲av | 久久久久久久精品吃奶| 久久精品成人免费网站| 国产国语露脸激情在线看| 欧美一级毛片孕妇| 在线观看免费高清a一片| 超碰成人久久| 每晚都被弄得嗷嗷叫到高潮| 国产主播在线观看一区二区| 视频区欧美日本亚洲| 日韩人妻精品一区2区三区| 成人黄色视频免费在线看| 狂野欧美激情性xxxx| 久久热在线av| 中文字幕av电影在线播放| 无限看片的www在线观看| 亚洲 欧美 日韩 在线 免费| 亚洲专区中文字幕在线| 啦啦啦在线免费观看视频4| 中文字幕精品免费在线观看视频| 老汉色∧v一级毛片| 亚洲成人精品中文字幕电影 | 黄色视频不卡| 老鸭窝网址在线观看| 最新在线观看一区二区三区| www.精华液| 中文字幕av电影在线播放| 又黄又粗又硬又大视频| 亚洲伊人色综图| 大陆偷拍与自拍| 久99久视频精品免费| 丁香六月欧美| 99久久人妻综合| 久久久久久久午夜电影 | 欧美黄色淫秽网站| 91大片在线观看| 日韩国内少妇激情av| a级毛片黄视频| 自线自在国产av| 岛国视频午夜一区免费看| 狠狠狠狠99中文字幕| 在线观看日韩欧美| 9色porny在线观看| 欧美成人午夜精品| 久久国产精品男人的天堂亚洲| 亚洲av日韩精品久久久久久密| 母亲3免费完整高清在线观看| 亚洲人成电影免费在线| 国产亚洲欧美精品永久| 免费在线观看日本一区| 超碰成人久久| 12—13女人毛片做爰片一| 国产成人系列免费观看| 在线观看免费视频网站a站| 欧美激情极品国产一区二区三区| 亚洲精品中文字幕一二三四区| 欧美精品一区二区免费开放| 最新在线观看一区二区三区| 在线观看免费高清a一片| 午夜91福利影院| 亚洲欧美激情在线| 最好的美女福利视频网| 亚洲成人国产一区在线观看| 国产精品自产拍在线观看55亚洲| 亚洲一码二码三码区别大吗| 亚洲一卡2卡3卡4卡5卡精品中文| 精品日产1卡2卡| 亚洲欧美日韩高清在线视频| 一本大道久久a久久精品| 热re99久久精品国产66热6| 欧美av亚洲av综合av国产av| 一进一出好大好爽视频| 国产精品自产拍在线观看55亚洲| 欧美乱妇无乱码| 国产在线精品亚洲第一网站| 美女国产高潮福利片在线看| 久久久国产欧美日韩av| 国产黄a三级三级三级人| 麻豆国产av国片精品| 国产高清videossex| 日韩精品中文字幕看吧| 男人舔女人下体高潮全视频| 亚洲国产精品一区二区三区在线| aaaaa片日本免费| 国产亚洲精品综合一区在线观看 | 两个人看的免费小视频| 国产精品免费视频内射| 亚洲成人精品中文字幕电影 | 久久天躁狠狠躁夜夜2o2o| 精品久久蜜臀av无| 91九色精品人成在线观看| 人人妻,人人澡人人爽秒播| 欧美日本中文国产一区发布| 青草久久国产| 国产人伦9x9x在线观看| 亚洲精品一卡2卡三卡4卡5卡| 精品国产一区二区三区四区第35| 亚洲性夜色夜夜综合| 国产欧美日韩一区二区精品| 精品国产乱子伦一区二区三区| 又紧又爽又黄一区二区| 国产xxxxx性猛交| 久久久久久大精品| 别揉我奶头~嗯~啊~动态视频| 国产在线精品亚洲第一网站| 在线观看免费视频日本深夜| 亚洲av美国av| 亚洲色图 男人天堂 中文字幕| 久久午夜综合久久蜜桃| 国产精品亚洲一级av第二区| 在线观看免费视频网站a站| 欧美激情高清一区二区三区| 欧美人与性动交α欧美精品济南到| 黑人猛操日本美女一级片| 欧美不卡视频在线免费观看 | 成年版毛片免费区| 国产成人系列免费观看| 亚洲人成77777在线视频| 在线观看日韩欧美| 制服人妻中文乱码| 一二三四在线观看免费中文在| 人人妻,人人澡人人爽秒播| 国产精品九九99| 两性夫妻黄色片| 午夜成年电影在线免费观看| 欧美 亚洲 国产 日韩一| 国产有黄有色有爽视频| a级毛片在线看网站| 欧美另类亚洲清纯唯美| 好看av亚洲va欧美ⅴa在| 亚洲av成人不卡在线观看播放网| 亚洲熟妇熟女久久| 麻豆国产av国片精品| 夜夜躁狠狠躁天天躁| 色哟哟哟哟哟哟| 天堂√8在线中文| 欧美日韩视频精品一区| 国产日韩一区二区三区精品不卡| 色哟哟哟哟哟哟| 日韩三级视频一区二区三区| 国产精品偷伦视频观看了| 久久国产精品人妻蜜桃| 国产亚洲av高清不卡| 免费高清视频大片| 国产国语露脸激情在线看| 91精品国产国语对白视频| 国产一区二区三区视频了| 欧美日韩中文字幕国产精品一区二区三区 | 日韩av在线大香蕉| 亚洲欧美精品综合久久99| 手机成人av网站| 精品国产美女av久久久久小说| 久久久久国内视频| 级片在线观看| 美女国产高潮福利片在线看| 免费看a级黄色片| 免费搜索国产男女视频| 精品第一国产精品| 黄色 视频免费看| 久久香蕉国产精品| 久久久精品国产亚洲av高清涩受| 一本大道久久a久久精品| 午夜福利一区二区在线看| 中文欧美无线码| 在线观看免费高清a一片| 久久久久亚洲av毛片大全| 免费搜索国产男女视频| 我的亚洲天堂| 亚洲av成人不卡在线观看播放网| 精品福利永久在线观看| 成人三级做爰电影| 亚洲精品一卡2卡三卡4卡5卡| 亚洲av成人av| 成年女人毛片免费观看观看9| 妹子高潮喷水视频| 亚洲性夜色夜夜综合| 黑人猛操日本美女一级片| 国产精品综合久久久久久久免费 | 亚洲五月色婷婷综合| 亚洲精品美女久久久久99蜜臀| 一进一出抽搐动态| 老司机亚洲免费影院| 一a级毛片在线观看| a级毛片在线看网站| 久久 成人 亚洲| 久久国产乱子伦精品免费另类| 亚洲黑人精品在线| 淫妇啪啪啪对白视频| 亚洲 国产 在线| 精品国产乱码久久久久久男人| 欧美中文综合在线视频| 妹子高潮喷水视频| 国产有黄有色有爽视频| 国产激情久久老熟女| 91字幕亚洲| 日韩av在线大香蕉| 看片在线看免费视频| 久久国产乱子伦精品免费另类| 欧美色视频一区免费| 欧美乱色亚洲激情| 国产精品久久久av美女十八| 亚洲国产精品一区二区三区在线| 91大片在线观看| 韩国精品一区二区三区| 午夜亚洲福利在线播放| xxx96com| 久久久久亚洲av毛片大全| 亚洲精华国产精华精| 国产激情欧美一区二区| 久久国产乱子伦精品免费另类| 麻豆av在线久日| 黄色毛片三级朝国网站| 国产真人三级小视频在线观看| tocl精华| 韩国av一区二区三区四区| 在线永久观看黄色视频| 日本一区二区免费在线视频| 成人永久免费在线观看视频| 婷婷六月久久综合丁香| 久9热在线精品视频| 最新在线观看一区二区三区| 国产精品久久久久久人妻精品电影| 窝窝影院91人妻| 精品乱码久久久久久99久播| 日本撒尿小便嘘嘘汇集6| 999久久久精品免费观看国产| 老鸭窝网址在线观看| 日日爽夜夜爽网站| 精品国产一区二区三区四区第35| 国产精品偷伦视频观看了| 韩国av一区二区三区四区| 国产主播在线观看一区二区| 无人区码免费观看不卡| 少妇 在线观看| 亚洲av电影在线进入| 亚洲男人的天堂狠狠| 老司机靠b影院| 免费人成视频x8x8入口观看| 人妻丰满熟妇av一区二区三区| 亚洲成人免费av在线播放| 亚洲人成电影观看| 亚洲avbb在线观看| 手机成人av网站| 亚洲九九香蕉| 黑人猛操日本美女一级片| 伊人久久大香线蕉亚洲五| 久久久久久久精品吃奶| 真人做人爱边吃奶动态| 丰满的人妻完整版| e午夜精品久久久久久久| 成人av一区二区三区在线看| 久久久久久久久久久久大奶| 久久精品亚洲精品国产色婷小说| 天堂动漫精品| 黑人猛操日本美女一级片| 久久久国产欧美日韩av| 亚洲一区二区三区欧美精品| 91精品国产国语对白视频| avwww免费| 18禁黄网站禁片午夜丰满| 9热在线视频观看99| 99热只有精品国产| 国产精品秋霞免费鲁丝片| 精品少妇一区二区三区视频日本电影| 亚洲中文av在线| 欧美另类亚洲清纯唯美| 一进一出抽搐gif免费好疼 | 亚洲精品国产色婷婷电影| 婷婷精品国产亚洲av在线| 中国美女看黄片| 搡老熟女国产l中国老女人| 亚洲熟妇熟女久久| 国产精品一区二区在线不卡| 国产精品久久久av美女十八| 国产熟女午夜一区二区三区| 亚洲av成人一区二区三| 欧美黑人欧美精品刺激| e午夜精品久久久久久久| 精品国产超薄肉色丝袜足j| 精品国产乱码久久久久久男人| 久久影院123| 在线视频色国产色| 麻豆国产av国片精品| 亚洲国产精品合色在线| 99精品欧美一区二区三区四区| 女人被狂操c到高潮| 黄色 视频免费看| 18禁裸乳无遮挡免费网站照片 | 精品高清国产在线一区| 巨乳人妻的诱惑在线观看| 国产成人一区二区三区免费视频网站| 久久久久久免费高清国产稀缺| 久久久久久久久久久久大奶| 久久精品亚洲熟妇少妇任你| 国产精品爽爽va在线观看网站 | 亚洲国产精品999在线| 免费在线观看黄色视频的| 亚洲自偷自拍图片 自拍| 亚洲中文av在线| 免费av中文字幕在线| 亚洲少妇的诱惑av| 后天国语完整版免费观看| 亚洲欧美精品综合一区二区三区| 亚洲狠狠婷婷综合久久图片| 国产精品一区二区在线不卡| 啦啦啦免费观看视频1| 亚洲人成伊人成综合网2020| 12—13女人毛片做爰片一| www.www免费av| 亚洲午夜精品一区,二区,三区| 亚洲性夜色夜夜综合| 国产1区2区3区精品| 88av欧美| 色尼玛亚洲综合影院| 搡老岳熟女国产| 亚洲欧美激情综合另类| 免费在线观看影片大全网站| 中亚洲国语对白在线视频| 老司机深夜福利视频在线观看| av有码第一页| 欧美久久黑人一区二区| 亚洲精品久久成人aⅴ小说| 国产精品亚洲av一区麻豆| 好男人电影高清在线观看| 如日韩欧美国产精品一区二区三区| 中文字幕精品免费在线观看视频| 国产免费现黄频在线看| xxxhd国产人妻xxx| 亚洲第一青青草原| 国产精品 欧美亚洲| 欧美在线黄色| 可以免费在线观看a视频的电影网站| 嫩草影视91久久| 欧美人与性动交α欧美精品济南到| 午夜两性在线视频| 亚洲国产欧美一区二区综合| 亚洲自偷自拍图片 自拍| 国产精品av久久久久免费| 一级毛片精品| 国产精品一区二区三区四区久久 | 手机成人av网站| 日本a在线网址| 精品久久久精品久久久| 国产成人精品在线电影| 国产熟女xx| 免费女性裸体啪啪无遮挡网站| 欧美乱色亚洲激情| 日韩欧美免费精品| 精品久久久精品久久久| 国产亚洲精品久久久久5区| 日韩免费高清中文字幕av| 天天添夜夜摸| 一本综合久久免费| 大陆偷拍与自拍| 美国免费a级毛片| 日日干狠狠操夜夜爽| 老熟妇乱子伦视频在线观看| 欧美日韩亚洲综合一区二区三区_| 亚洲av片天天在线观看| 日韩视频一区二区在线观看| 男男h啪啪无遮挡| 视频区欧美日本亚洲| 欧美最黄视频在线播放免费 | 男人的好看免费观看在线视频 | 欧美日韩亚洲高清精品| 欧美色视频一区免费| 国产一区二区三区视频了| 波多野结衣一区麻豆| 国产精品综合久久久久久久免费 | 在线免费观看的www视频| 天堂影院成人在线观看| 久久中文看片网| 最新美女视频免费是黄的| 亚洲成人国产一区在线观看| 黑人欧美特级aaaaaa片| 十八禁网站免费在线| 一级毛片精品| 视频区欧美日本亚洲| 精品国产国语对白av| 99国产综合亚洲精品| 又黄又粗又硬又大视频| 69av精品久久久久久| 久久香蕉精品热| 国产精品偷伦视频观看了| 亚洲国产欧美日韩在线播放| 别揉我奶头~嗯~啊~动态视频| 国产精品影院久久| 大码成人一级视频| 美女大奶头视频| 欧洲精品卡2卡3卡4卡5卡区| 99国产精品免费福利视频| 色综合婷婷激情| 久久久久久久精品吃奶| 日韩中文字幕欧美一区二区| 麻豆国产av国片精品| 成人国产一区最新在线观看| 日韩一卡2卡3卡4卡2021年| 国产三级在线视频| 自拍欧美九色日韩亚洲蝌蚪91|