PAGE 1
i SYNCHRONIZATION OF C OMPLEX COUPLED OSCIL LA TORS UTILIZING CONSENSUS CONTROL by NATHAN DANIEL CLARK B.S., University of Colorado Denver, 2013 A thesis submitted to the Faculty of the Graduate School of the University of Colorado in partial fulfillment of the requirements for the degree of Master of Science Electrical Engineering 2014
PAGE 2
ii This thesis for the Master of Science degree by Nathan Daniel Clark has been approved for the Electrical Engineering Program by Miloje Radenkovic Chair Mohsen Tadi Yimeng Deng November 13, 2014
PAGE 3
iii Clark, Nathan (M.S., Electrical Engineering) Synchronization of Complex Coupled Oscillators Utilizing Consensus Control Thesis directed by Professor Miloje Radenkovic ABSTRACT Synchronization phenomena has found its roots in many different scientific and engineering disciplines that have s ought to understand the phenomena of how oscillating systems coupled through some medium are able to synchronize in their steady state. Modeling synchronization of coupled oscillators is highly complex and nonlinear and is an ongoing area of rich research and development. Synchronization has been established in such circles as biology, medicine, signal processing, dynamical system control and many other areas. This thesis contributes to further the understanding of how to maintain control of inherently synchronous systems that may be perturbe d by outside disturbanc es through development of an algorithm that has previously never been utilized in this context This algorithm is applied to complex coupled oscillators as proposed by Yoshiki Kuramoto. Today these complex coupled oscillators, also known as Kuramoto Oscillators, are used to describe any system of nearly identical osc illators that are tenuously coupled. Another assumption that Kuramoto oscillators makes is that each oscillator has its own natural frequency associated with it and any exchange of information between two oscillators are sinusoidally related by the phase difference between them.
PAGE 4
iv This thesis concludes with a brief application of consensus contr ol to a weakly connected system of networked power system elements, as well as, a broad survey of the numerous applications that are currently being explored. In addition, this thesis provides a first fundamentals background and introduction to Kuramotos oscillators and consensus control. The form and content of this abstract are approved I recommend its publication. Approved: Dr. Miloje Radenkovic
PAGE 5
v DEDICATION I ded icate this thesis to my wife, To nya Waldschmidt, who without her support and constant encouragement, I would never have been able to complete my education. I want to thank her for sacrificing all the time we missed together over the past 6 years which allowed me the opportunity to be fully engaged in my education and perform to the best of my ability. Tr uly, without her constant love and support I would not be where I am today.
PAGE 6
vi ACKNOWLEDGMENT I would like to sincerely extend my gratitude to my advisor Dr. Miloje Radenkovic, for whom his support and guidance have been invaluable to my entire career aspirations and educational development. His commitment to mentoring and teaching me over the past year has been a major contributing factor to my engineering capacities, thorough understanding of complex subject matter and support to seeing this thesis to its final form I want to commend Dr. Radenkovic for his hard work and time commitment to both always striving to be an exceptional teacher of engineering while also continuously driving the el ectrical enginee ring program towards excellence
PAGE 7
vii CONTENTS Chapter 1. Introduction ................................................................................................................. 1 1.1 Purpose and Objective ................................................................................................ 6 1.2 Thesis Layout ............................................................................................................... 6 2. Overvie w of Kuramoto Oscillators ............................................................................. 8 2.1 Importance of the Critical Coupling Constant ............................................................. 9 2.2 Deriving the Coupling Constants Relationship to Modulus Order Parameter r ..... 11 2.3 Determination of the Critical Coupling Constant ...................................................... 12 2.4 Further Expansion on Kuramoto Oscillators ............................................................. 16 2.5 Another Physics Based Complex Network Metric ..................................................... 18 3. Consensus Control ........................................................................................................ 21 3.1 Graph Theory and Networks ..................................................................................... 22 3.2 Discrete time Average Consensus Algorithm ............................................................ 24 4. Fitting Consensus Control to Complex Oscillator Networks ..................................... 27 4.1 Adaptive Control Preliminaries: A Normalized Gradient Algorithm ......................... 28 4.2 Normalized LMS B ased Consensus Algorithm ........................................................... 31 4.3 Simulation Experiment .............................................................................................. 35 5. Conclusion ................................................................................................................. 47 Bibliography ...................................................................................................................... 49 Appendix A. Kuramotos Algorithm MATLAB Program .................................................................. 51 B. Non Linear Consensus Control .................................................................................. 54 C. Consensus Control of Power Ne twork Modeled by Kuramoto Oscillators ............... 56 D. Consensus Control of Po wer Network Matpower Exchange .................................... 58
PAGE 8
viii LIST OF FIGURES F igu r e 2 1: Order parameter denoting the cohesiveness of a particular system ................. 11 2 2: Phase relationships of 10 Kuramoto oscillators ...................................................... 13 2 3: Frequency relationships of 10 Kuramoto oscillators ............................................... 14 2 4: Order parameter as a function of coupling strength for 10 oscillators ................... 15 2 5: Phase relationhsips of 50 Kuramoto oscillators ...................................................... 16 2 6: Frequency relationships of 50 Kuramoto oscillators .............................................. 17 2 7: Order paramete r as a function of coupling strength for 50 oscillators ................... 17 3 1: Example of how a Laplacian graph is constructed from its degree and adjacency matrix ................................................................................................................................ 23 3 2: Example topo logies for different example linkages for consensus control ............ 26 4 1: Process flow diagram detailing the interaction of consensus and adaptive control with Kuramoto oscillators ................................................................................................. 33 4 2: Simulation of 10 oscillators with random initial frequencies .................................. 34 4 3: Process flow diagram detailing the information exchange between MatPower and the consensus algorithm ................................................................................................... 39 4 4: Cumulative phase angles for IEEE 6 bus power network system ........................... 40 4 5: System parameter estimates with an abrupt change in power demanded ............ 41 4 6: Frequency and power demanded/generated for the IEEE 6 bus system ................ 42 4 7: Perturbed generator frequencies and power load generator responses ............... 43 4 8: Cumulative phase angles of a 57 bus power network system ............................... 44 4 9: Parameter estimates of a 57 bus power network system ....................................... 45 4 10: Frequencies of the 57 bus power network syste m at varying times .................... 46
PAGE 9
1 1. I ntroduction The fundamental concept of synchronization can be traced as far back as the 17th century Dutch mathematician Christian Huygens. Huygens had observed that two clocks with pendulums hanging on a wall wou ld eventually synchronize in motion. Upon closer inspection Huygens realized that the primary reason that these two clocks would synchr onize was due to the fact that the imperceptible motion of each of the pendulums was causing a motion in the beam upon which they were attached. The beam would cause a feedback motion on each of the pendulums until they would be mutually synchronized and 180 degrees out of phase with each other. At this point, the dual feedback from each of the pendulums that were now 180 degrees out of phase but completely frequency synchronized caused a netting effect of forces applied to the beam and the beam would cease to move anymore. At this point the system of the beam and the two clocks were in a steady state that Huygens referred to as the symphony of two clocks [1]. Huygens description of synchronization through the symphony of two clocks provides a deep int uitive understanding to the phenomenon of synchronization among two or more coupled oscillators. Through the weak coupling of the wooden beam the two clocks would transfer their kinetic energy from the motion of the pendulums to the beam, and the beam wou ld consequently transfer its relative motion to each of the two clocks and their respective pendulums. This discovery of mutual synchronization striving for harmony and steady state has found interest in many disciplines today such as biology and nature, e cology, sociology, engineering, neuroscience and arts.
PAGE 10
2 Scientists within the field of biology took notice to the phenomenon of synchronization when the y noticed the effect of the sun on the rising and falling of leaves throughout a 24 hour period. Synchronization phenomenon may not seem all that apparent in an observation of an organism that is moving itself for radiation in the process of photosynthesis to sustain itself, but it is intriguing when the scientists place the same plants in a dark room and als o notice the same rising and falling motion of the leaves. Their concl usion was that the plants were still synchronizing to other indirect effects which the sun was placing on their surrounding environment [1]. Additionally, biologists hav e analyzed the phenomenon of synchronization on the molecular level with genetic networks. One of the most studied examples is that of a synthetic gene network called a repressilator. The repressilator is a network of three genes that repress each other during transcription in a cyclic way [16]. The concentration levels through time of each of these three genes can be modeled by, [ ] = [ ] + 1 + [ ] = ( [ ] [ ]) where [] is the mRNA concentration encoded by [] is the conce ntration of the translated protein is the promotor rate and is the ratio of the protein decay rate to the mRNA decay rate. This system has a unique steady state and has drawn the interest of biologists in understanding biological rhythms amongs t thousands of cellular oscillators coordinating in synchronism.
PAGE 11
3 Ecologists have also uncovered a natural 10 year cycle in lynx populations within three regions of Canada. These scientists discovered that the lynx populations were in synchrony with the mig ration patterns that they exhibited in finding food. These foodwebs as there known in the field of ecology are pervasive across nature whenever there is interplay between vegetation, herbivores and predators [16]. Synchronization is also quite prevalent in the study of neural systems where there are billions of neurons connected together in a network. It has been found that complex neuronal networks play a significant role in cognitive processing. In fact, when the brain is processing information it does so in the most economic possible way. The brain minimizes the global cost of metabolism and energy to create new connections. Therefore most connections and synapses are local in nature although there are special ins tances where there are longer connections necessary. By and large the brain will attempt to synchronize with local neurons to the greatest extent possible [16]. Another example of synchronization can be found in the synchronization of flashing fireflies. T he consideration that flashing fireflies also follow this same phenomenon is an important one because it suggests that any formulaic model of synchronization can be applied on a large scale to many oscillators in this case, fireflies. Fireflies when obse rved flash their organic lights in synchrony and as other fireflies join the larger group they also join in unison. This self organization can also be observed with applauding audiences most commonly noted in Eastern Europe [2] At first audiences applaud out of synch and the clapping
PAGE 12
4 sounds distorted and noisy. After several oscillations of clapping, the applause tends to become synchronized after a select few lead individuals produce enough consistent applaud in unison continuously. The lead indiv iduals are necessary for the followers to join in on the synchronized clapping. This concept of leader follower will be discussed in more detail later on in this thesis. S ynchronization can be observed and leveraged in computer science and engineering disciplines. In data mining, which is the process of searching for patterns in large amounts of data to extract meaning and information, computer scientists encode data vectors representative of a database as frequencies of oscillators. The hope is that s ynchronization of the different frequencies will group similar data together and therefore categorically give meaning to a heterogeneous space of complex data [16]. Computer scientists have proposed the following dynamical model, ( ) = ( ) + sin ( ) ( ) where is a multivariate data point, ( ) is the l th component of the phase vector, is a function that determines the proximity of the phase vector to the distance = between each data point, is a c oupling constant and is the number of data points. Finally, s ynchronization can even be found in the most unlikely of places such as the dispersal of opinion formation among humans in social sciences. The idea that synchronization could drive opinion formation was born from the desire to understand
PAGE 13
5 how long it would take for partial or complete consensus to emerge out of a set of different opinions. One of the most i nteresting models proposed takes into account the fact that two individuals with very strong differing opinions will tend to interact less with each othe r. In other words, their coupling strength decreases. This is modeled by exponentially tenuous opinions of a particular individual to another individuals opinions. Through time each of these individuals will cease communication because they eac h perceive the ongoing relationship with diminishing returns [16]. Social scientists have pro posed the following model, = + sin where is the intrinsic opinion of an individual, and and are scaling factors that determine the coupling strength. R ather the cultural glue that biases certain opinions into agreement or not. It is clear that there are a plethora of examples of synchronization to be found throughout the literature. This brie f survey is to provide a wide appreciation of how prevalent synchronization can be used to model any number of phenomena that exhibits a state transition from heterogenic complex systems to those of homogenous systems be it through spatial and/or temporal transitions.
PAGE 14
6 1.1 Purpose and Objective The purpose and obje ctive of this thesis is to build upon the understanding of synchronization and how it relates to control system techniques such as consensus control. Beyond simply relating the two schools of thought it is the intent of this thesis to marry the two approaches so that they complement each other. Kuramoto Oscillators are employed to model the system dynamics of numerous complex oscillators that are initially out of synchronization both in phas e and frequency. The thesis then lays the framework for how to apply the technique of consensus control to maintain synchronization of complex oscillator networks. The stability of the oscillators communicating through the Kuramoto model and controlled through consensus is then stress tested or more specifically Perturb ed and Observe d to evaluate the effectiveness and robustness of the Kuramoto model in conjunction with the control algorithm. 1.2 Thesis Layout Chapter 2 gives a brief overview of Yoshiki Kuramot os classic coupled oscillators. Some insight is provided into preliminaries necessary to appropriately use this model as well as a discussion about the strength of coupling necessary between each oscillator in order for the oscillators to converge to a co mmon frequency. Chapter 3 introduces the concept of consensus control and gives historical examples of its use in practice. In addition, the consensus control algorithm is put into the context of
PAGE 15
7 large networked oscillators that are weakly and strongly co nnected and communicate directionally and bidirectionally. Chapter 4 delves into the framework and methods proposed by this thesis to effectively integrate consensus control into Kuramoto Oscillators for the purpose of controlling weakly coupled systems and/or easily bifurcated systems. The successful implementa tion of the algorithm on control is proved through perturb and observe methods that seek to drive a synchronized power system out of both frequency and phase synchrony. Chapter 5 concludes the thes is with final observations and thoughts towards future research and tangible imple mentation of our new algorithm.
PAGE 16
8 2. Overview of Kuramoto Oscillators Yoshiki Kuramoto first proposed his model of complex oscillators and how best to model their behavior towards in phase and out of phase synchrony in his 1975 contribution to Lecture Notes in Physics, International Symposium on Mathematical Problems in Theoretical Physics [3] Kuramoto was not the first to produce a mathematical model representing coupled oscillators. Arthur Winfree is actually attributed to laying down the framework for modeling the behavior of large interacting oscillators. Although his model is not often used his assumptions and insight into the correct implementation of the model still stand today. Winfree proposed that for the oscillators to achieve synchrony they need to first be nearly identical and second that their coupling is weak. These two assum ptions were necessary to reduce the complexity of the math describing a large network of complex coupled oscillators that would otherwise have many state variables. Winfrees assumptions allowed the model to be reduced to describing each oscillators effec t on the group at large through their phase relationships to each other. Kuramoto took Winfrees assumptions and further expanded on them by allowing Winfrees model to be generalized among any type of network ed system of oscillators. Kuramoto also provide d weighting factors to the phase relationship between each oscillator to better model the bounds necessary for synchrony. It was these contributions that enabled Kuramotos oscillator model to gain academic notoriety [2]. Kuramotos equation is given as:
PAGE 17
9 = + asin = 1 , (2.1) Or in the special case of = 1 = + sin = 1 , (2.2 ) where ( ) are the coupling parameters is the total number of networked oscillators and is the oscillator connected to sinusoidal phase relationship of oscillators and the oscillators natural frequency For simplicity in modeling it is often best practice to set the initial frequencies of normally distributed about a center controlling mean frequency of the modelers choice. As mentioned previously, the main draw from researchers to the Kuramoto model is in the constant or rather more commonly known as the critical coupling strength [4]. Kuramoto proved that if could exceed a particular threshold constant in relation to the distribution of the natural fr equencies than synchronization would be certain. 2.1 Importance of the Critical Coupling Constant Kuramoto determined that there is value for which a particular coupling constant is < When this is the case it is not possible for a system of networked oscillators to achieve synchrony. However, when the coupling constant exceeds ( > ) the critical coupling constant than individual oscillators are able to achieve synchron ization
PAGE 18
10 A very useful way to visually describe Kuramotos model is through set of nodes on a unit circle distributed randomly around the circle initially and given random initial frequencies. Kuramoto modeled these nodes on a unit circle to help understa nd the critical coupling constant. This new performance measure he introduced, he coined as the performance measure of phase cohesiveness known as the complex order parameter. Kuramoto defined this order parameter as: = (2.3 ) where is the modulus of the order operator or measure of synchronization and is argument of the order. Another way to view the order operator is to say that it represents the mean collective behavior of all the networked coupled oscillators. In other words, when = 0 the system is said to be completely unsynchronized, whereas when = 1 the system is said to be in complete harmony and synchrony. This order parameter can be seen as an arrow pointing out towards the edge of the circle with its magnitude represented by the synchrony of the system as shown in figure 2 1
PAGE 19
11 Figure 21 : Order parameter denoting the cohesiveness of a particular system As can be seen from figure 2 1 the length of r increases as the oscillators converge together in both phase and frequency synchronism. 2.2 Deriving the Coupling Constants Relationship to the Modulus Order Parameter r It is beneficial to relate the coupling constant to the modulus order parameter to gain a better insight into how the modulus parameter, which is a metric of synchrony, relates to the strength of coupling We begin by multiplying both sides of equation 2. 3 by ( )= = 1 , (2.4 ) considering only the imaginary parts we are left with sin ( ) = sin = 1 , (2.5 )
PAGE 20
12 Equation 2.4 can now be substituted into equation 2.1 to arrive at a form that represents the Kuramoto model in terms of the modulus parameter and the coupling constant: = + ( ) = 1 , (2.6 ) Equation 2.5 offers us useful insight into the synchronization of complex oscillators that was not obvious from equation 2.1. As we can see, the oscillators synchronize towards a mean order parameter It is also clear that with small coupling constants of that each oscillator will adhere to its own natural frequency On the other hand with a large value of each oscillator feels the influence of its neighbors either pushing or pulling and will entrain together toward the mean order paramete r [4]. 2.3 Determination o f the Critical Coupling Constant Determining the order parameter as it relates to the coupling constant provides us with intrinsic dynamics about a particular system. When we observe the order parameter, which is the meas ure of synchrony, against a varying coupling constant it allows us to identify a knee, or rather the critical coupling constant that must be greater than to guarantee synchronism. To portray this we consider a system with 10 oscillators distributed randomly both in their initial phase and frequency. Figure 2.2 illustrates the phase relationship of the 10 oscillators through time with a varying coupling constant of 0 < < 3 5
PAGE 21
13 Figure 2 2 : Phase relationships of 10 Kuramoto oscillators As figure 2 2 shows there is relatively good phase synchronization of 8 of the 10 oscillators. The two oscillators that have drifted away from the group have synchronized in frequency as shown as in figure 2 3 0 20 40 60 80 100 120 140 160 180 200 14 12 10 8 6 4 2 0 2 4 TimeTheta
PAGE 22
14 Figure 2 3 : Frequency relationships of 10 Kuramoto oscillators What the previous two figures illustrate is that phase synchronism and frequency synchronism are independent of each other. While the 10 oscillators did achieve frequency synchronism they did not achieve phase synchronism. As noted earlier, for there to be complete phase and frequency synchronism the modulus order parameter must be equal to one. In reality can never reach one but rather approaches one as as is suggested by figure 2 4. 0 10 20 30 40 50 60 70 5 4 3 2 1 0 1 2 3 4 TimeFrequency
PAGE 23
15 Figure 2 4 : Order parameter as a function of coupling strength for 10 oscillators With only 10 Kuramoto oscillators it has become evident that very good simultaneous phase and frequency synchronization may not be possible depending on the initial phase and frequenc y of each of the oscillators. However, the power of this analysis is evident in that we can now determine that the oscillators tend to be more synchronized at coupling constants of 3 or greater. 0 0.5 1 1.5 2 2.5 3 3.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Coupling strength: KAverage Coherence: r
PAGE 24
16 2.4 Further Expansion on Kuramoto Oscillators As shown in figure 2 2 and 2 3 we can see that the drifting oscillators are almost mirror images of each other and therefore cancel each other out in their contribution to the modulus order parameter From statistics we know that outliers are eventually phase d out as we increase the number of samples taken from a population. This phenomenon is well known as the Central Limit Theorem. In much the same way, Kuromoto Oscillators exhibit this same type of behavior as more identical oscillators are introduced to t he complex networked system. Figure 2 5, 2 6 and 2 7 clearly demonstrates how increasing the amount of oscillators in a network has the effect of improving dynamic non linear behavior by increasing the tightness of each of the oscillators phase and frequ ency relationships to the group as a whole. Figure 2 5 : Phase relationhsips of 50 Kuramoto oscillators 0 20 40 60 80 100 120 140 160 180 200 3 2.5 2 1.5 1 0.5 0 0.5 1 1.5 2 TimeTheta
PAGE 25
17 Figure 2 6 : Frequency relationships of 50 Kuramoto oscillators Figure 27 : Order parameter as a function of coupling strength for 50 oscillators 0 20 40 60 80 100 120 4 3 2 1 0 1 2 3 TimeFrequency 0 1 2 3 4 5 6 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Coupling strength: KAverage Coherence: r
PAGE 26
18 The critical coupling constant for this particular system is still three but it is clear that the effect of many more oscillators added to a particular system has the effect of smoothing the coupling strength contribution to the average coherence of the system, as well as tightening the phase relationships between all the oscillators. 2.5 Another Physics Based Complex Network Metric The critical coupling strength metric is useful for determining when a complex network will synchronize because the critical coupling constant is a measure of pulling and pushing strength between each of connected nodes within a network. This metric, however, is not very useful for evaluating the connectedness of a particular complex network. One useful metric for determining the connectedness or betweeness of a network is a metric u sed most commonly in the physics community known as the clustering coefficient. The clustering coefficient calculates the number of shortest paths through any given node to help understand how communicative the adjacency matrix is [16]. This clustering coe fficient formulation is known as, = = ( ) (2.6) where is the number of nearest neighbors connected to and is the pr obability degree distribution of node A large clustering coefficient implies that there are many paths within the network and many redundant paths to get between any two nodes. A low clustering coefficient is analogous to a sparse graph is which there are not many paths and therefore the network is not very well connected.
PAGE 27
19 We have described in detail the analytical approach as given by Kuramoto to model synchronism of complex oscillators. Kuramoto s model is, however, intrinsically a continuous time model with information exchange between each of the nodes as modeled as smooth sinusoidal phase couplings. There are instances where Kuramotos model is not very useful for example where the coupling between nodes and their neighbors is through pulses. T hese types of models are known as pulse coupled models [16]. Pulse coupled models are used to represent nodes within a network that have linear phase dynamics. However, once a certain node h as reached a threshold that node then sends a pulse (delta fun ction) to the other nodes within the network and subsequently relaxes back to a state of rest. Pulse coupled models may be better understood through the scaling relations between the time it takes to achieve synchronization and the topological parameters o f a network. The relationship takes the form as, ~ where is the time needed to achieve synchronization, is the number of links, is the number of nodes and and are constants generally set at 1.3 and 1.5 respectively. This mod el shows that the time needed to synchronize a random network inc reases with randomness and sparsity Pulse coupled models are worthy of mention because synchronization is seemingly always associated with continuous time systems. Pulse coupled models are in fact a rich
PAGE 28
20 and intriguing way to model how the brain fires off neurons and how different synapses a re made within the brain when learning and muscle movement activities are engaged.
PAGE 29
21 3. Consensus Control Consensus control is understood to be a network of agents, or nodes, within a particular system that strives to reach agreement among all the connected agents to a certain state of interest [5]. Consensus control is a special form of an optimization problem that seeks to find a steady state solution with out relying on an objective function. Rather consensus control seeks harmony and agreement among all connected agents through communication to other agents. With that said, it should be obvious that there can be direct corollaries made between the phenomenon of Kuramoto oscillators, consensus control and synchronization. The refore it is a natural fit to suppose that Kuramotos model and the method of consensus control are complimentary to each other. In fact, as will be discussed in chapter 4, the framework will be laid out that relates consensus control to Kuramotos model. Consensus co ntrol is a heuristic problem solving tool that takes observations of nature and tries to fit it to a mathematical model. Examples of a heuristic mathematical model is that of particle swarm optimization that was brought up in the scientific community to be tter understand and model how large groups of birds flocking together can move in unison and not collide into each other. From these observations of flocking birds it was noted that all the birds did not perfectly move in unison when their observations were recorded and slowed down. Rather, there was a delay of propagation from a global leader bird to the last bird to react to what eventually the flock was doing together.
PAGE 30
22 Scientists who noted this delay realized that each bird was reacti ng to what its near est neighbors where doing. I n turn all the other birds would react the same way and so on and so forth. With that, it was relatively easy to come up with a simplified mathematical model to mimic what th e synchronized flocking birds w ere doing [15 ] 3.1 Graph Theory and Networks G raph theory is a useful tool to describe the topology of a network Each element within the graph/matrix = ( ) represents an agent/node/element = { 1 2 , } that contributes to the whole of the system. Each of the nodes are connected by thei r edges denoted as All edges that connect to a particular agent/node and its neighbors is given as = ( ) = [ ] The graph representing all the interconnected agents can be summarized as follows, ( ) = ( ) ( ) + ( ) ( 0 ) = (3.1) To better understand the topology of a particular graph is helpful to define a Laplacian graph which is a graphs degree matrix minus its adjacency matrix. The degree matrix denotes the number of neighbors connected to a graph along the diagonal whereas the adjacency matrix describes how each agent is physically connected. An example of what the Laplacian graph might lo ok like is shown in figure 3 1
PAGE 31
23 Figure 3 1 : Example of how a Laplacian graph is constructed from its degree and adjacency matrix With this understanding the dynamics of the agents can be modeled as = = [ ] (3.2) The agents of the L aplacian graph are given as = 1   = (3.3) The value of consensus is intuitive for network of agents in a graph. It is simply ( ) = 1 / where as mentioned previously corresponds to the initial state of each of the agents We can now see that for a well connected graph, that is to say that we can get to any other node through all the edges that there is a mean value that incorporates the initial state of each of the agents [15]. Some important rema rks are in order now that pertain to this thesis in particular. The first of which is that of the concept of cooperation. For all agents to be in cooperation each individual agent must be willing to participate and allow its objective to update to the grou ps overall global objective. One of the main goals of this thesis is not only to have cooperation among many different agents but to have control over the global objective. To do this we employ the leader follower method that selects a random agent
PAGE 32
24 to be a leader and all other agents are followers. This notion is in direct contradiction to the requirement that all agents be cooperative and allow their objective to be flexible. Introducing a leader that is not willing to change its objective actually does n ot compromise the global objective of the networked system. It is only when there are two or more leaders with opposing objectives that this becomes a separate problem and issue and gives rise to the rich subject of game theory which is beyond the purpose and scope of this thesis. Within the literature there is also mention of the concept of the synchronization by pacemaker which is in spirit what the leader follower method is. The pacemaker is usually a single/sub set of nodes that do not change their frequency and therefore entrain the rest of the network to follow and synchronize. The pacemaker functions only within a network where the phases of each node are organized in a hierarchical fashion. This is to say that nodes are grouped together at nearly identical phases at different distances from the pacemaker within the network [16]. 3.2 Discretetime Average Consensus Algorithm The discrete time average consensus algorithm is a slight variation from the previously discussed general consensus algorithm. This algorithm functions primarily based on a notion of particles moving through a virtual space with specific headings and direc tions. Each particle updates its heading and speed based upon the average of its own heading and velocity and that of its nearest neighbors. Much like the Kuramoto oscillators, which
PAGE 33
25 is also a variation on consensus, the group of particles will eventually find agreement on direction and heading [6]. For agents in an undirected network topology we can say that there exists an asymptotically stable average consensus that follows the form, ( + 1 ) = ( ) ( ) (3.4) where = the number of a gents present within the system and ( ) = [ ( ) , ] is a vector that relates the agent to all of its neighbors. We can say that to have equal weights of applied across all agents allows us to propose the following, ( ) = (3.5) Equation 3.4 can now be rewritten in the following form, ( + 1 ) = ( ) (3.6) which assumes that each agent can influence any other agent present within the networked system even if a particular agent does not have many neighbors. The graph must be connected and directed. When we say the graph must be connected, mathematically speaking we are saying that the matrix is irreducible or rather that the rank is 1 A good way to visualize co nnectedness is through figure 3 2
PAGE 34
26 Figure 3 2 : Example of topologies of different example linkages for consensus control In figure 3 2 the first and second topologies are not strongly connected T here is not a directed path to every agent. With the first topology we can only send information from agent 3, to agent 2, to agent 1. Agent 3 does not exactly know what agent 1 is doing. Likewise for the second topology, agent 1 cannot fully know what is happening with agent 3. However, the last topology clearly illustrates that all agents can be sought out. Agent 1 may still not exactly k no w what agent 3 was doing in the previous time step but for the purposes of consensus control it is no t necessary since agent 3 will eventually be informed by its previous information even though its original information has been aggregated and maybe distorted through agents 1 and 2. Consensus control has been briefly described in its most simplistic form to allow the reader to get a feel for its intrinsic and intriguing characteristics. Consensus control is widely understood and utilized in industry in such areas as; military UAV drones, distributed sensor networks and air traffic control to name a few. The end goal of introducing consensus control is to give a brief introduction of the specific averaging algorithm to Kuramoto oscillators which will be described in detail in the follow ing c hapter.
PAGE 35
27 4. Fitting Consensus Control to Complex Oscillator Networks It has been suggested several times throughout this thesis that there exists a natural corollary between Kuramotos oscillators and the consensus algorithm. We also note that these two well studied models can be used in tandem to craft a much more robust control algorithm that can be applied to a number of different models as introduced in chapter 1 It is the purpose of this chapter to lay the framework and process to successfully im plement a novel control algorithm that utilizes Kuramotos oscillator model and consensus control from which is the main contribution of this thesis. Synchronization of mutually coupled oscillators is achieved through several key criteria: 1) Each of the oscillators must be identical in their behavior or nearly identical. 2) Each of the oscillators must exhibit a natural frequency. 3) Each of the oscillators must be coupled through some sort of medium that seeks a steady state. 4) There exists a critical coupling constant relating each of the different oscillators upon when achieved and exceeded guarantees synchronization. The four aforementioned points are what enable synchronization to take place. In reality we know that rarely are networked systems in such a state of homogony with the exception of harmonious systems within nature. Manmade systems on the other hand such as telecommunicati ons and power systems are always in a critical balance of state between synchronous convergence and inadvertent times uninte ntionally disturbed
PAGE 36
28 divergence. This critical balancing point for both linear and non linear systems is referred to as a bifurcation point in a mathematical context. Bifurcation is essentially when a system at rest is perturbed into two separate and diver gent states. In dynamical studies small changes in a particular parameter could have the effect of completely transitioning a particular system from its original topological state to a completely different one. Since systems that are naturally described through Kuramotos model are always at the brink of unrest and bifurcation it makes sense to wrap a control algorithm around this model to both give more stability to the complex oscillators and to allow the system to quickly and efficiently re organize i tself when disturbed. 4.1 Adaptive Control Preliminaries : A Normalize d Gradient Algorithm Before we can delve into the main algorithm that this thesis is based upon we need to spend some time developing a simple control algorithm that will be used to estimate the system parameters. Without parameter estimation there is no way for the c onsensus algorithm to recover the Kuramoto oscillators from some exter nal disturbance or even know what state the system was at to begin with. A gradient adaptive technique is discussed, that through time the system parameters are always known and the consensus algorithm can focus on bringing in all agents to a common objective. Another way to look at it is to say that without knowledge of the system and space each agent cant po ssible now where they are going much like telling a large
PAGE 37
29 group of people to find a single exit in a dark room. Those people are going to bump into each other a nd not be successful in finding the exit and eventually get frustrated with each other and diverge and regress in frustration. To begin we let a system be described by the fo llowing difference equation 4.1 ( + 1 ) + ( ) + + ( + 1 ) = ( ) + ( 1 ) + + ( + 1 ) where ( ) is the system input and ( ) is the system output. We can represent ( 1 ) unit delay with the operator ( ) = ( 1 ) which gives us, ( 1 + + + ) ( + 1 ) = ( + + + ) ( ) (4.2 ) Each side of 4.2 can be also represented in shorthand notation if we allow ( ) = 1 + + + and ( ) = + + + ( ) where ( ) / ( ) is referred to as the transfer operator. Furthermore we can define a set of sy stem parameters as = [ , ; , ] and a signal vector as, ( )= [ ( ) , ( + 1 ) ; ( ) , ( + 1 ) ] 4.1 can then be written as, ( + 1 ) = ( ) (4.3 )
PAGE 38
30 We assume that the parameters are unknown and make it our main goal to determine these parameters from observations of t he signal and output vectors to arrive at an estimate of which is denoted as Our final realization of our system for adaptive control is now given as, ( + 1 ) = ( + 1 ) ( ) (4.4 ) With our system described we can now extrapolate a gradient algorithm that will be used to assist in synchronizing Kuramoto oscillators by adaptively finding a particular system s parameter estimates. This part is important sinc e we cant control the system if we dont know the parameters describing the system which can change if the system is disturbed. With that being said, we can estimate the parameters of a system by taking the gradient of the error of the system with respect to The error of the system is given as, ( )= ( ) ( 1 ) (4.5 ) taking the gradient of this we get equation 4.6 ( )= ( 1 ) ( ) ( 1 ) (4.6 ) The direction of this gradient is one in which a small change in will yield a an appreciable change in ( ) With this consideration, we are inc lined to make the recursive algorithm as follows, ( ) = ( 1 ) + ( ) ( 1 ) ( ) ( 1 ) ( 1 ) (4.7 )
PAGE 39
31 where is a scalar that contributes to the rate of convergence of the parameter estimator. We can further refine by normalizing it with respect to the square of the signal vector This is given as, ( ) =  ( )  (4.8 ) where 0 < < 2 and the number one in the denominator is to assure that there is no division by zero. The algorithm is known as the Kaczmans or projection algorithm. Inserting 4.8 into 4.7 we have our final recursive adaptive algorithm. ( ) = ( 1 ) +  ( )  ( 1 ) ( ) ( 1 ) ( 1 ) (4.9 ) 4.2 Normalized LMS Based Consensus Algorithm With the conclusion of the adaptive control method presented we are now able to construct a novel algorithm and framework to reliably control and maintain synchronism among complex oscillators. To begin we know that we can relate the an oscillators discrete frequencies to phases as follows, ( + 1 ) =( ) ( ) (4.10) ( + 1 ) = ( ) + ( + 1 ) (4.11) Where are the frequencies, are the phases and the freque ncies and phases are related to each other by ( + 1 ) = with being the time step Noting that
PAGE 40
32 Kuramotos model is a relationship between a node s frequency and the summation of all its connected nodes related sinusoidally through their phase s we can formulate equation 4.11 as, ( ) ( ) = ( ) ( 0 ) = (4.12) w here ( ) = ( + 1 ) sin ( + 1 ) ( + 1 ) ( ) sin ( ) ( ) (4.13 ) a nd where are the known couplings between each oscillator. Further we can represent 4.12 as shown in 4.14 which is a form better suited for parameter estimation as, ( + 1 ) = ( ) + ( ) ( ) (4.14) where = and ( ) = ( ) To control and update ( ) we employ the average consensus algorithm along with the leaderfollower method as, ( ) = ( ) (4.15 ) where is unique in that it describes only the set of a gents currently connected to the node. To maintain knowledge of the space we need to adjust with the adaptive normalize d gradient algorithm as d erived earlier and repeated here for convenience as, ( + 1 ) = ( ) ( ) ( ) [ ( + 1 ) ( ) ] (4.16)
PAGE 41
33 w here initial ( 0 ) are arbitrary parameter estimates. Computing each of these equations fo r a number of time steps described through the known couplings will allow ( 0 ) agents to s ynchronize. Figure 4 1 illustrates this process through a high level flow diagram. Figure 4 1 : Process flow diagram detailing the interaction of consensus and adaptive control with Kuramoto oscillators
PAGE 42
34 The aforementioned algorithm was derived based on the work by Radenkovic and Tadis Self tuning Synchronization Algorithm [18] Figure 4 2 illustrates the algorithm functioning on a general sys tem of oscillators with randomly d istributed initial frequen cies and phases about zero and with a single leader hol ding its frequency constant at zero Figure 4 2 : Simulation of 10 oscillators with random initial frequencies It is clear from figure 4 2 that only within 11 steps through the algorithm that almost complete synchronization is achieved and that 9 of the 10 oscillators are entraining towards the leader held at zero. 1 2 3 4 5 6 7 8 9 10 11 1 0.5 0 0.5 1 1.5 2 2.5 Time [steps]Frequency
PAGE 43
35 4.3 Simulation Experiment As mentioned previously the application to be simulated as Kuramoto Oscillators is that of a power network system. The goal of this section of the thesis is to prove that in a simulated environment that the algorithm models reality as close as possible and proves reliable and autonomous control of a po wer system with con sensus control and gradient algorithms. We selected to apply the algorithm to a small power system with a few buses, also known as a microgrid, to prove that the algorithm works. A microgrid provides a good test bed for the algorithm because it exhibits certain characteristics that would cause it to fall out of synchrony. Microgrids are essentially a collection of generators and loads that are in perfect balance and connected through physical means of buses, also known as nodes, and transmission lines known in graph theory context as edges. Since electricity generation and consumption is cyclical in nature and requires that all generation to be totally consumed, it is a natural fit to be modeled by Kuramotos oscillators and controlled through consensus and a simple gradient algorithm. We begin by defining our network model and fit it to Kuramotos oscillator model. It is well known that a generator and load can be modeled through its dynamic physical equations known as swing equations. The swing equations presented here are given as is and are not derived since it is beyond the purpose and scope of this thesis. The well known generator and load swing equations are given as, = , { 1 , } (4.1 7 )
PAGE 44
36 where there are generators/loads, each modeled by their inertia constant electrical power output mechanical power input d amping constant rotor/phase angle and the synchronously rotating frame With the nodes defined we still need to know the topology of the network. For this we define the network through its edge/adjacency matrix represented physically by an a dmittance matrix which describes both the physical connections of each of the nodes/agents of the system and resistances through them. We can then use the fact that an admittance matrix represents all the connected conductances and susceptances, which are merely the real and imaginary part of resistance in transmission lines. With these facts we have that the electrical power output as, = cos + sin (4.18) The and are the conductance and susceptance between the and generator with internal voltages Substituting 4.18 into 4.17 we have, = , sin (4.19 ) where is the feeback gain. This second order nonlinear model is far too complex and unwieldy to analyze. We therefore put it into an approximate form, also known as the non uniform Kuramoto model. The non uniform Kuramoto model serves as the bridge between 4.19 and Kuromotos oscillator model as given by 2.1 To derive the non uniform Kuramoto model we first define the natural frequency of the agents as
PAGE 45
37 = which is the mechanical power input into a generator, the coupling between each oscillator as = which is the maximum power transferred between each generator, a perturbation parameter = / and the speed of convergence represented as = ( / ) /( /) From these definitions we can represent 4.19 as, = + sin (4.20) In conventional power systems is rather small in 60 Hz environments. Equation 4.20 then reduces to, = sin (4.21) With further algebraic cancellation and choosing the feedback gain to be nearly zero for further simplification, we find the final form of the modeled power system network in its most nearest form to a Kuramoto Oscillator as, = sin (4.22) With insight established between the relationship of a power network system has to Kuramoto s model we can begin to draw a corollaries between Kuramotos non uniform model to a network of generators and loads The first system simulated is an IEEE 6 bus power network system that models a small system connected to a utility that has two generators and three loads. The loads are
PAGE 46
38 modeled after the active power drawn from the system that results in a non uniform Kuramoto power balance equation given as [4] + = sin ( ) (4. 23 ) where = =     ( ) for the specific case as it rela tes to a power system and as the inverse of which is modeling the damping coefficients of the power network system. We are going to assume that all the generators are connected to DC/AC inverters and that the deviation from its nominal power output > 0 is directly related to the change in frequency This results in the droop controlled inverter dynamics [11] = sin ( ) (4. 24) What is most striking about these two equations is how similar they are. Essentially the load and generator equations should balance each other out. In other words, what is inje cted into the system from 4.24 is what should be consumed by 4.23 Power flow in a complex power system network is rather complex and beyond the scope of this thesis. We therefore leveraged Matpower and free power flow program developed by Cornell University to execute the proper power flow equations for us within Matlab. Figure 4 3 gives a process flow diagram to detail how the consensus algorithm exchanges information to and from Matpower.
PAGE 47
39 Figure 4 3 : Process flow diagram detailing the information exchange between MatPower and the consensus algorithm
PAGE 48
40 With the framework established between Matpower and our algorithm we tested the IEEE 6 bus system against it to obs erve its performance. Figure 4 4 demonstrates the phase angles are in consens us and not deviating Figure 4 4 : Cumulative phase angles for IEEE 6 bus power network system Figure 4 5 shows how the adaptive normalized LMS algorithm immediately updates when there is a change in the power demanded by the loads at time step 80. As can be seen the algorithm is rather good at readjusting itself to accommodate the new system dynamics. 0 20 40 60 80 100 120 10 0 10 20 30 40 50 60 70 Time [steps]Cumulative Phase Angles Utility Gen1 Gen2 Load1 Load2 Load3
PAGE 49
41 Figure 4 5 : System parameter estimates with an abrupt change in power demanded Figure 4 6 gives a side by side comparis on of how the algorithm is able to maintain synchrony even after there has been a rapid change in power demand. The algorithm first brings all agents to the global objective of 60 Hz and then maintains synchrony until there is an abrupt change in power dem anded. As illustrated there is a cool down period for the algorithm to recalculate the new system parameters and find the lead bus again within the system. 0 20 40 60 80 100 120 140 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Time [steps]Theta Hat System Parameter Estimates Utility Gen1 Gen2 Load1 Load2 Load3
PAGE 50
42 Fi gure 4 6 : Frequency and power demanded/generated for the IEEE 6 bus system Perturbing the loads gives insight into how the algorithm can self regulate and maintain synchronism among all agents within the system. This algorithm maintains synchrony even if the generators themselves are pertur bed as illustrated in figur e 4 7 0 20 40 60 80 100 120 55 60 65 Time [steps]Frequency [Hz] Utility Gen1 Gen2 Load1 Load2 Load3 0 20 40 60 80 100 120 500 0 500 Time [steps]Power [W] Utility Gen1 Gen2 Load1 Load2 Load3
PAGE 51
43 Figure 4 7 : Perturbed generator frequencies and power load generator responses Here we can see that the frequencies of the generators are perturbed at time steps 20, 40, and 80 and we can see how that affects the loads and generators power curves. We can also extend the control algorithm to include additional buses. We selected a 57 bus system to test the robustness of the algorithm under much more complex co nditions. Figure 4 8 is the phase relationship of all 57 buses. 0 20 40 60 80 100 120 50 55 60 65 70 Time [steps]Frequency [Hz] 0 20 40 60 80 100 120 500 0 500 Time [steps]Power [W] Utility Gen1 Gen2 Load1 Load2 Load3 Utility Gen1 Gen2 Load1 Load2 Load3
PAGE 52
44 Figure 4 8 : Cumulative phase angles of a 57 bus power network system Figure 4 9 provides insight into how the parameter estimates are updated as the system is perturbed at time steps 40 and 80. 0 20 40 60 80 100 120 140 10 0 10 20 30 40 50 60 70 80 Time [steps]Cumulative Phase Angles
PAGE 53
45 Figure 4 9 : Parameter estimates of a 57 bus power netwo rk system Finally we show the effect of perturbing the 57 bus system when we perturb slightly at time step 5 and perturb greatly at time steps 40 and 80. Figure 4 10 illustrates how the algorithm is still able to maintain frequency at 60 Hz. 0 20 40 60 80 100 120 140 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Time [steps]Theta Hat System Parameter Estimates
PAGE 54
46 Figure 410 : Frequencies of the 57 bus power network system at varying times 0 20 40 60 80 100 54 56 58 60 62 64 Time [steps]Frequency [Hz]
PAGE 55
47 5 Conclusion This thesis main goal was to introduce Kuramoto oscillators in their raw form as Yoshiko Kuramoto modeled them originally in 1 975, as well as, introduce the concept of consensus control. Through a preliminary understanding of Kuramotos oscillators we were able to show that for a certain constant coupling between all interconnected oscillators there is a value that which all oscillators will synchronize this is known as the critical coupling constant. The critical coupling constant was explored to show that even though Kuramotos oscillators model the phenomenon of natural synchronization in nature, not all oscillators in reality are able to achieve harmony. Even if modeled oscillators in reality are able to achieve synchronism, there exist external forces that could potentially disturb the system and cause the network oscillators to diverge and lose synchronism due to the fact that the system may be near its bifurcation point. In addition, this thesis explored the notion that consensus control wrapped around Kuramotos oscillators could provide both additional reliable synchron ism on Kuramotos model but with the added control through the leader follower concept in whic h one agent is maintained at a global objective. Finally, we were able to ascertain that there exists a very close variation of Kuramotos model that can be extended to a power system network with loads and inverters. This is important as it extends the theory to a potential application that the presented algorithm could be applied to. Through numerous examples of perturbing the power
PAGE 56
48 network system, we observed that the system is still able to maintain s ynchronism within a period of time. The main contribution of this thesis is in crafting the framework and implementation of a new algorithm that adds an additional layer of control around Kuramotos celebrated model through consensus control and an adaptiv e gradient algorithm.
PAGE 57
49 BIBLIOGRAPHY [1] A, Pikovsky, M. Rosenblum, J. Kurths, Synchronization: A Universal Concept in Nonlinear Sciences. Cambridge University Press: Nonlinear Science Series 12, 2001 [2] S.H. Strogatz, From Kuramoto to Crawford: Exploring the onset of synchronization in populations of coupled oscillators. Physica, 2000 [3] Y. Kuramoto, H. Araki, Lecture Note in Physics, International Symposium on Mat hematical Problems in Theoretical Physics. Springer, New York, 1975 [4] F. Dorfler, F. Bullo, Synchronization in Complex Oscillator Networks: A Survey. Automatica, 2014 [5] R. Olfati Saber, A. Fax, R. Murray, Consensus and Cooperation in Netwo rked Multi Agent Systems. Proceedings of the IEEE, Vol. 95, No. 1, 2007 [6] A. Jadbabaie, J. Lin, A.S. Morse, Coordination of Groups of Mobile Autonomous Agents Using Nearest Neighbor Rules. Automatic Control IEEE, Vol. 48, Issue 6, 2003 [7] N. Ainsworth, S. Grivalva, A Structure Preserving Model and Suffficient Condition for Frequency Synchronization of Lossless Droop Inverter Based AC Networks. IEEE Transactions on Power Systems Vol. 28, Issue 4, 2013 [8] W. Ren, R.W. Beard, Consensus s eeking in multi agent systems under dynamically changing interaction topologies. IEEE Trans. On Automatic Control, 2004 [9] J. Lin, A.S. Morse, B.D.O Anderson, The multiagent rendezvous problem. Proc. Of CDC, 2003 [10] H.G. Tanner, G.J. Pappas, V. Kumar, leader to formation stability. IEEE Trans. Robot. Autom., Vol. 20, No. 3, 2004 [11] J. Porco, F. Dorfler, F. Bullo, Synchronization and Power Sharing for Droop Controlled Inverters in Islanded Microgrids, Automatica, 2013 [12] S.H. Strogatz, Sync. Hyperion New York, 2001 [13] P. Ashwin, J.W. Swift. The dynamics of n weakly coupled identical oscillators. Journal of N onlinear Science, 1992 [ 1 4] B. C. Daniels, Synchronization of Globally Coupled Nonlinear Oscillators: The Rich Behavior of the Kuramoto Model Stanford University, 2005
PAGE 58
50 [15] J. Wang, D. Wang, Particle Swarm Optimization with a Leader and Followers. Progess in Natural Science, Elsevier, 2008 [ 16] A. Arenas, A. Diaz Guilera, J. Kurths, Y. Moreno, C. Zhou, Synchronization in Complex Networks. Physics Reports, Elesevier, 2008 [17] S.H. Strogatz, Exploring Complex Networks, Nature 410, 2001 [18] M. Radenkovic, M. Tadi, Self tuning Average Consensus in Complex Networks submitted for publication, 2013.
PAGE 59
51 A. Kuramotos Algorithm MATLAB Program % Simulation of Kuramoto's model (Chapter 6, Box A) % Written by Joakim Munkhammar % clear; % N is i, the number of individuals in % the collective (N>1 required)) N=50; % Total number of timesteps T=200 % The range for coupling constant K % Setting up the vector of Ktot Klow = 0; Kup = 6; Kstep = 0.05; totsteps = round((KupKlow)/Kstep); Ktot(1) = Klow; for i=2:totsteps; K tot(i) = Ktot(i1)+Kstep; end %Time step size. As tau > 0 approximation of continumm model improves tau=0.1 % w is defined as w(i), time independent % Initial condition for w, initial distribution % Standard is normal distribution (Scaled) but % Centered about zero. Unimodal with infinite tail. % Initial conditions for initial distribution % The distribution parameter gamma gamma=1; % The width of the distribution lower=10; upper=10; for run=1:10 run % Initial conditions with normal distribution for w sigN=1; w = random('Normal',0,sigN,1,N);% Random on initial conditions % The bifurcation loop for Steps=1:totsteps
PAGE 60
52 %Steps; % Setting the value of K each lap K=Klow+ Steps*Kstep; % theta is defined as theta(time,i) % Here initial conditions for theta for i=1:N theta(1,i)=w(i); end % Attempting to calculate r (Phase coherence between 0 and 1) % and phi (Average angle) for the first steps of time rx=0; ry=0; phi(1)=0; for i=1:N phi(1) = phi(1) + (1/N)*theta(1,i); % Calculation mean angle phi rx=rx+(1/N)*cos(theta(1,i)); % Sum of mean xpart of theta ry=ry+(1/N)*sin(theta(1,i)); % Sum of mean ypart of theta end r(1) = sqrt(rx*rx + ry*ry); r(2) = r(1); phi(2)=phi(1); % The main loop for t=2:T % Initial conditions each timestep in the loop rx=0; ry=0; phi(t+1)=0; % The loop of individuals for i=1:N % Main equation theta(t,i) = theta(t1,i) + tau*(w(i) + K*r(t)*sin(phi(t)theta(t1,i))); rx=rx+(1/N)*cos(theta(t,i)); % Sum of mean xpart of theta ry=ry+(1/N)*sin(theta(t,i)); % Sum of mean ypart of theta % Calculating mean angle phi for next step of time phi(t+1) = phi(t+1) + (1/N)*theta(t,i); end r(t+1) = sqrt(rx*rx + ry*ry); % Calculating total mean radius end % Assigning the final value of r at each K to rtot rtot(run,Steps) = mean(r(T100:T+1)); storew(Steps,:)=w + K*r(t)*sin(phi(t)theta(t1,:)); end end %pdf: g(0) and g''(0)for Normal distribution g0=1/(sigN*sqrt(2*pi)) % Plotting the radius r (coherence) over K as bifurcation diagram
PAGE 61
53 figure(1) hold on plot(Ktot,mean(rtot),'k')%b hlx=xlabel('Coupling strength: K'); hly=ylabel('Average Coherence: r'); axis([Klow Kup 0.05 1]) figure(2) %% Plotting inititial distribution of w(i) subplot(2,1,1) hist(w,[4:0.05:4]) axis([4 4 0 50]) title('Kuramotos model'); hlx=xlabel('i'); hly=ylabel('g(w) probability density'); %% Plotting final distribution of w(i) subplot(2,1,2) hist(storew(46,:),[4:0.05:4]) axis([4 4 0 50]) hlx=xlabel('i'); hly=ylabel('g(w) probability density'); figure (3) plot(storew) x label('Time') ylabel('Frequency') figure (4) plot(theta) xlabel('Time') ylabel('Theta')
PAGE 62
54 B. Non Linear Consensus Control clc clear all %Make Aij matrix A = [0 0 1 0 0 0; 0 0 0 0 1 0; 1 0 0 1 0 0; 0 0 1 0 1 0; 0 1 0 0 0 1; 1 0 0 0 1 0]; %Number of Particles N_i = 6; % A = ones(N_i,N_i); % A = [0 1 1 1 1 1; % 1 0 1 1 1 1; % 1 1 0 1 1 1; % 1 1 1 0 1 1; % 1 1 1 1 0 1; % 1 1 1 1 1 0]; %Gain Constant mu = 1; %Period Constant T_s = 0.01; %Initial position of the particles delta = 1*randn(1,N_i)'; %Initial parameter estimates theta_hat = [1 0 0 0 0 0]'; %Initial x (frequency) x = [59.5 59.8 60 60.1 60.2 60.7; 59.5 59.8 60 60.1 60.2 60.7]'; steps = 400; for t = 1:steps for i = 1:N_i delta(i,t+1) = delta(i,t) + T_s x(i,t); for j = 1:N_i delta_mem(i,j) = A(i,j)*(sin(delta(j,t+1)delta(i,t+1)) sin(delta(j,t)delta(i,t)))*T_s; end x(1,t) = 60; phi(i,t) = sum(delta_mem(i,:)); x_bar_avg = 0; for j = 1:N_i x_bar_avg = x_bar_avg + A(i,j)*x(j,t); end x_bar_avg = x_bar_avg + x(i,t); x_bar(i,t) = (1 / (1 + nnz(A(i,:)))) x_bar_avg; phi(i,t) = (x_bar(i,t) x(i,t)) (1 + nnz(A(i,:))); x(i,t+1) = x(i,t) + theta_hat(i,t) phi(i,t); theta_hat(i,t+1) = theta_hat(i,t) mu (phi(i,t) (x(i,t+1) x_bar(i,t))) / (1 + phi(i,t)^2); end end figure(1)
PAGE 63
55 plot(delta') legend('1','2','3','4','5','6') figure(2) plot(theta_hat') figure(3) plot(x')
PAGE 64
56 C. Consensus Control of Power Network Modeled by Kuramoto Oscillators clc clear all %Make Aij matrix A = [1 0 1 0 1 0; 0 1 1 0 1 0; 1 0 1 1 0 0; 1 0 1 1 1 0; 0 1 0 0 1 1; 1 0 1 0 1 1]; %Number of Particles N_i = 6; A = ones(N_i,N_i); % A = [0 1 0 1 1 0; % 0 0 1 1 1 1; % 0 0 0 0 1 1; % 0 0 0 0 1 0; % 0 0 0 0 0 1; % 0 0 0 0 0 0]; %Gain Constant mu = 0.5; %Period Constant T_s = 1; %Initial position of the particles delta = randn(1,N_i)'; %Initial parameter estimates theta_hat = [1 1 1 1 1 1]'; %Initial x (frequency) x = [59.9 59.9 60 60.1 60.2 60.7; 59.5 59.8 60 60.1 60.2 60.7]'; %x = [60 60 60 60 60 60; 60 60 60 60 60 60]'; freqs = x(:,1); steps = 300; for t = 1:steps for i = 1:N_i delta(i,t+1) = delta(i,t) + T_s x(i,t); for j = 1:N_i delta_mem(i,j) = A(i,j)*(sin(delta(j,t+1)delta(i,t+1)) sin(delta(j,t)delta(i,t)))*T_s; end phi(i,t) = sum(delta_mem(i,:)); x_bar_avg = 0; for j = 1:N_i x_bar_avg = x_bar_avg + A(i,j)*x(j,t); end x(6,t) = 60; x_bar_avg = x_bar_avg + x(i,t); x_bar(i,t) = (1 / (1 + nnz(A(i,:)))) x_bar_avg; phi(i,t) = (x_bar(i,t) x(i,t)) (1 + nnz(A(i,:)));
PAGE 65
57 x(i,t+1) = x(i,t) + theta_hat(i,t) phi(i,t); theta_hat(i,t+1) = theta_hat(i,t) mu (phi(i,t) (x(i,t+1) x_bar(i,t))) / (1 + phi(i,t)^2); end if t == 1 opt = mpoption; opt = mpoption(opt,'OUT_ALL',0,'VERBOSE',0); mpc = runpf('case6ww',opt); [power_shared, freqs, A, mpc1, delts] = final_power_sharing (x(:,t),A,t,mpc,delta,theta_hat,T_s); P_i(:,t) = [power_shared(1:3) power_shared(4:6)]; else [power_shared, freqs, A, mpc1, delts] = final_power_sharing (x(:,t),A,t,mpc1,delta,theta_hat,T_s); P_i(:,t) = [power_shared(1:3) power_shared(4:6)]; end if var(x) > 0.000001 ff = 0; else x(:,t+1) = freqs; end delta(:,t+1) = delts; if t == 40 x(2,t+1) = 59.97; end if t == 60 x(2,t+1) = 60.97; end if t == 80 x(3,t+1) = 60.06; end end figure(1) plot(delta') legend('Utility','Gen1','Gen2','Load1','Load2','Load3') figure(2) plot(theta_hat') legend('Utility','Gen1','Gen2','Load1','Load2','Load3') %axis([0 50 1 1]) figure(3) subplot(2,1,1) %x = x/norm(x); plot(x','LineWidth',1.5) grid legend('Utility','Gen1','Gen2','Load1','Load2','Load3') axis([0 300 50 70]) subplot(2,1,2) plot(P_i','LineWidth',1.5) grid legend('Utility','Gen1','Gen2','Load1','Load2','Load3') axis([0 300 1000 1000])
PAGE 66
58 D. Consensus Control of Power Network Matpower Exchange Note: It is necessary to have Appendix D and C sitting in the same folder as Matpower for the algorithm to function correctly. function [power_shared, x, A,mpc1, delta] = final_power_sharing ( freqs, Aij, t, mpc, delta, theta_hat, T_s ) opt = mpoption; opt = mpoption(opt,'OUT_ALL',0,'VERBOSE',0); %Convert Frequency in Hz to Radians/s freqs = 2*pi*freqs; %Extract useful data from MatPower define_constants bus_type = mpc.bus(:,BUS_TYPE); PowD = mpc.bus(:,PD); PowDind = find(PowD); for ii = 1:length(PowDind) PD(ii,1) = PowD(PowDind(ii,1),1); end PG = mpc.gen(:,PG); %Make Aij for branch flows from_bus = mpc.branch(:,F_BUS); to_bus = mpc.branch(:,T_BUS); from_bus_power = mpc.branch(:,PF); to_bus_power = mpc.branch(:,PT); Aij_From_Flows = full(sparse(from_bus,to_bus,from_bus_power,6,6)); Aij_To_Flows = full(sparse(to_bus,from_bus,to_bus_power,6,6)); Aij_Branch_Flows = Aij_From_Flows Aij_To_Flows; %Get deltas at each bus delta = pi/180*mpc.bus(:,VA); %Define the damping coefficient as inverse of theta_hats D_i = 1./theta_hat; %Define the Invertor and Load Kuramoto Swing Equations for power for i = 1:6 for j = 1:6 PFr(i,j) = Aij_Branch_Flows(i,j)*sin(delta(j)delta(i)); end P_from(i) = sum(PFr(i,:)); end P_i(1) = D_i(1)*freqs(1) + P_from(1); P_i(2) = D_i(2)*freqs(2) + P_from(2); P_i(3) = D_i(3)*freqs(3) + P_from(3); P_i(4) = 1 (D_i(4)*freqs(4) + P_from(4)); P_i(5) = 1 (D_i(5)*freqs(5) + P_from(5)); P_i(6) = 1 (D_i(6)*freqs(6) + P_from(6));
PAGE 67
59 %Input Updated Bus Data Back into MatPower mpc.gen(1,2) = abs(P_i(1)); mpc.gen(2,2) = abs(P_i(2)); mpc.gen(3,2) = abs(P_i(3)); mpc.bus(4,3) = 70; mpc.bus(5,3) = 70; mpc.bus(6,3) = 70; if t < 200 mpc.bus(4,3) = 0.1*abs(P_i(4)); mpc.bus(5,3) = 0.1*abs(P_i(5)); mpc.bus(6,3) = 0.1*abs(P_i(6)); else mpc.bus(4,3) = 0.51*abs(P_i(4)); mpc.bus(5,3) = 0.35*abs(P_i(5)); mpc.bus(6,3) = 0.1*abs(P_i(6)); end %Rerun power flow in Matpower and extract new power flows, angles and frequencies mpc = runpf(mpc,opt); define_constants bus_type = mpc.bus(:,BUS_TYPE); PowD = mpc.bus(:,PD); PowDind = find(PowD); for ii = 1:length(PowDind) PD(ii,1) = PowD(PowDind(ii,1),1); end PG = mpc.gen(:,PG); delta = pi/180*mpc.bus(:,VA); %Make Aij for branch flows from_bus = mpc.branch(:,F_BUS); to_bus = mpc.branch(:,T_BUS); from_bus_power = mpc.branch(:,PF); to_bus_power = mpc.branch(:,PT); Aij_From_Flows = full(sparse(from_bus,to_bus,from_bus_power,6,6)); Aij_To_Flows = full(sparse(to_bus,from_bus,to_bus_power,6,6)); Aij_Branch_Flows = Aij_From_Flows Aij_To_Flows; %Define the Invertor and Load Kuramoto Swing Equations for frequency for i = 1:6 for j = 1:6 PFr(i,j) = Aij_Branch_Flows(i,j)*sin(delta(j)delta(i)); end P_from(i) = sum(PFr(i,:)); end freqs(1) = abs((P_i(1)P_from(1))/D_i(1)); freqs(2) = abs((P_i(2)P_from(2))/D_i(2)); freqs(3) = abs((P_i(3)P_from(3))/D_i(3)); freqs(4) = abs(1*(P_i(4)+P_from(4))/D_i(4)); freqs(5) = abs(1*(P_i(5)+P_from(5))/D_i(5)); freqs(6) = abs(1*(P_i(6)+P_from(6))/D_i(6));
PAGE 68
60 %Return Frequency, Adjacency Matrix and Delta back to consensus algorithm power_shared(1) = mpc.gen(1,2); power_shared(2) = mpc.gen(2,2); power_shared(3) = mpc.gen(3,2); power_shared(4) = mpc.bus(4,3); power_shared(5) = mpc.bus(5,3); power_shared(6) = mpc.bus(6,3); mpc1 = mpc; x = freqs/(2*pi); A = zeros(6,6); for ii = 1:6 for jj = 1:6 if Aij_From_Flows(ii,jj) ~= 0 A(ii,jj) = 1; end end end
