Iza Marfisi

Games in Higher Education

Keywords: Game-based Learning, Serious Game, Learning Games, Simulation Games

Introduction

This entry presents an overview of how and why Learning Games are used in higher education.

Learning Games can be defined as games that are designed to captivate the learners' attention and facilitate their learning process. They have explicit educational purposes and can be used for teaching at all levels of education. All types of games can be used for learning: board games, card games, roleplaying games, First Person Shooter games, simulation games, management games, puzzle games, treasure hunts…

The main characteristic of Learning Games for higher education is the fact that they are designed to teach specific complex skills taught at university or during professional training programs.

Unfortunately, it is not infrequent to observe strong opposition on the part of this target audience to this mode of learning, that these adult students associate with children.

The use of Learning Games in primary school seems natural to teachers and is encouraged by specialists in didactics and neuroscience. This learning technique is much less frequently used in middle school and is almost completely absent from higher education. Yet teachers at all these levels are faced with the same problems, such as lack of motivation and investment, for which games are known to be an effective solution. This entry presents an overview of the games that can be used for higher education and the reasons why some teachers and students still show resistance to this type of learning. The numerous advantages of games for higher education will then be presented, citing games presently used in universities, in graduate schools and for professional training. Finally, this Draft : Marfisi-Schottman I. (2019) Games in Higher Education. In: Tatnall A. (eds) Encyclopedia of Education and Information Technologies. Springer, Cham entry presents the current research questions that need to be addressed concerning the design of games for higher education and the acceptance of these games by teachers.

Accepting Games in Higher Education

Teaching with games

As children, we naturally use games as a pedagogical tool to enhance our emotional, sensory, motor, cognitive, intellectual and social development. This natural process is adopted by teachers in primary schools to teach mathematics, history, art, music and foreign languages [START_REF] Kamii | Group Games in Early Education: Implications of Piaget's Theory[END_REF]. This pedagogical approach aims to use game mechanics in order to captivate students' attention and engage them in their own learning process [START_REF] Dondlinger | Educational video game design: A review of the literature[END_REF]. When playing, the learners become the central actors of their learning process, a role quite unlike the passive position they occupy, most of the time, in traditional education. In order to win the Learning Game, students need to make decisions based on their newly acquired knowledge. Moreover, the learners find themselves emotionally engaged in the game and this facilitates the memorization of their decisions (National Research Council, 2000). Games can therefore be used to facilitate learning certain skills, for which traditional teaching methods are not satisfactory [START_REF] Mayo | Games for science and engineering education[END_REF].

What games for teaching?

These past years, many digital Learning Games have been developed. Digital Learning Games do not require any physical material and offer many advantages such as adapting to the learner's profile, allowing access to the game anytime, anywhere and displaying infinite patience when it comes to repeating concepts. In addition, the new generation of students are accustomed to playing computer games. An American study carried out in 2000 showed that, on average, students who have obtained an undergraduate degree spent only 5,000 hours reading compared to 10,000 hours playing video games and 20,000 hours watching television [START_REF] Prensky | Digital Natives, Digital Immigrants[END_REF]. With the democratization of smartphones and tablets over the last 20 years, this tendency has no doubt been amplified. In 2015, 75% of Americans between 13 and 17 had access to smartphones and more than half had access to tablets [START_REF] Lenhart | Teens, Social Media & Technology Overview[END_REF]. In this context, it is in our interest to provide tools and methods to help teachers use the attractiveness of video games to facilitate the acquisition of complex skills.

In order to be effective, Learning Games need to offer engaging game mechanics that serve the educational objectives. There are many theories on the best way to choose and integrate these game mechanics with the learning content. The extrinsic motivation approach, for example, was very common during the 80s and 90s, when Learning Games were referred to as "edutainment". With this approach, inspired by the behaviorist learning theories, gaming and learning are treated as separate entities. Indeed, the learners are usually first asked to perform certain educational exercises and, if they succeed, they are given access to a short game or receive points and badges as a reward. This approach met much criticism, and was often regarded as "sugar coating" over drill and practice exercises [START_REF] Bruckman | Can Educational Be Fun?[END_REF]. Nevertheless, this approach has become quite popular in higher education these last years under the term "gamification". Indeed, it is easy to gamify existing courses on online platforms used by universities and MOOCs (Massive Open Online Courses) [START_REF] Gené | Gamification in MOOC: Challenges, Opportunities and Proposals for Advancing MOOC Model[END_REF]. Many of these platforms offer plugins to create badges, progress bars, leader boards and they show fun content only if the learners have achieved a certain goal (a certain number of points or activities…). In regards to the limitations of this approach, several Learning Game designers have chosen to explore another direction called intrinsic motivation approach [START_REF] Fabricatore | Learning and Videogames: An Unexploited Synergy[END_REF][START_REF] Kafai | Game Design as an Interactive Learning Environment for Fostering Students' and Teachers' Mathematical Inquiry[END_REF], influenced by the constructivist learning theories. It consists in blending game mechanics with the learning content instead of using them separately. The idea is to choose a game design that is adapted to the educational goal and weave it into the pedagogical activities. According to numerous studies [START_REF] Habgood | Zombie Division: Evaluating principles of intrinsic integration[END_REF][START_REF] Lepper | Intrinsic motivation and instructional effectiveness in computerbased education, Aptitude, Learning and Instruction[END_REF][START_REF] Ryan | Self-determination theory and the facilitation of intrinsic motivation, social development, and well-being[END_REF], it is this cohesion between learning and fun that truly engages the students in their activity and facilitates the learning process. However, this approach implies creating custom Learning Games in collaboration with teachers and game designers, a very costly process that is not always possible. By combining both approaches, maximum efficiency can be obtained.

It is important to understand that the effectiveness of a Learning Game depends not only on the characteristics of the Learning Game artifact, but also on the way it is used [START_REF] Sanchez | A praxeological perspective for the design and implementation of a digital role-play game[END_REF].

Therefore, the game mechanics and the educational content do not necessarily need to be embedded in the Learning Game artefact itself. For example, teachers can use basic exercises in a fun context by distributing rewards, setting up battles between groups or introducing the concept with a story… The other extreme is also possible: using a pure game, such as Assassin's creed1 , for educational purposes (teaching History2 ).

Resistance to Learning Games

The use of Learning Games in the context of higher education is subject to many forms of resistance.

First of all, the use of games with adult or young adults can be problematic because they often view games as futile and only fit for children [START_REF] Eyster | Playing Game in Class, A beginner's guide for Higher Education[END_REF]. This opinion is sometimes shared by colleagues and superiors. Teachers who use Learning Games therefore need to prepare a flawless argumentation that shows how their Learning Games will help students reach the given educational goals.

Another critical concern is the student's evaluation: is it fair to evaluate learners according to their scores and actions in the Learning Game? The nature of games, as defined by [START_REF] Caillois | Man, Play and Games[END_REF] implies that playing can only be done free willingly and the actions taken in the game should not have any consequences on the real world. Using a Learning Game to grade students therefore contradicts the very notion of play.

Finally, it is unrealistic to believe that one Learning Game will suit all students. We all have different player profiles (Seeker, Survivor, Daredevil, Mastermind, Conqueror, Socializer and Achiever)3 and therefore do not necessarily like the same types of games [START_REF] Nacke | BrainHex: A neurobiological gamer typology survey[END_REF]. This must be kept in mind when designing Learning Games: it is best to choose game mechanics that cover several player profiles.

Even though the use of Learning Games faces resistance in higher education and requires complex multidisciplinary design skills, many have proven their efficiency in various domains. In the next section, the advantages that games can offer for higher education will be presented and illustrated with examples of games presently used in universities, graduate schools and professional training programs.

Advantages of Games for Higher Education

An opportunity to simulate and manipulate

Computer games have the advantage of offering adaptable virtual environments that are very useful for recreating specific situations and simulating the context in which learners will use their skills. This is particularly advantageous when the context is impossible or very difficult to reproduce because of its costly or dangerous nature. For example, The Resuscitation Game4 is used to teach medical procedures and reanimation technics for neonatal resuscitation. Rail Simulator5 (Figure 1) allows future train drivers to validate their management and driving skills. Virtual Reality is one of the latest advances in simulations. This technology immerses the learners in an interactive virtual word and allows them to practice their orientation and manipulation skills. Virtual reality is also used for simulating stressful situations in order to help learners control their emotions [START_REF] Marfisi-Schottman | Enabling Teachers to Create Authentic Interview Simulations[END_REF][START_REF] Ponder | Immersive VR decision training: telling interactive stories featuring advanced virtual human simulation technologies[END_REF].

Enhanced emotional engagement

Learning Games have the power of involving students in their own learning process with game mechanics such as competition, rewards, social recognition -and many more -that enhance motivation and activate the students' capacities [START_REF] Dondlinger | Educational video game design: A review of the literature[END_REF]. Learning Games offer much more than a simple simulation environment. They offer the possibility of creating role playing games based on complex scenarios with stories and quests. Thanks to these mechanics, learners can project themselves into a character and are emotionally engaged in the action of helping this character attain his/her goal. Throughout the game, the character will need to master various skills to fulfill his/her objectives. Many educational games offer such engaging scenarios. Starbank The Game6 (Figure 2) for example is used to train new bank employees, Les aventures de Casey Warren7 is used to train employees in data security.

Relief from routine

Using games is a good way to break the routine that sets in after several months of class. Introducing this teaching method changes the usual dynamics of the class and gives a second chance to the students who feel they are too far behind or that have given up all interest in the class. The games also require a set of new skills (communication, organization, tactical planning) that are not usually put forward in class. This provides teachers with fresh insights into their students as learners [START_REF] Kirriemuir | Literature Review in Games and Learning[END_REF]. Games also create a break, a state of "relaxed attention" that enhances learning [START_REF] Thiagarajan | Design your own games and activities: Thiagi's templates for performance improvement[END_REF]. CheckiO 8 (Figure 3), for example, offers a new approach for learning Python and JavaScript programming languages.

A new student-teacher relationship

Games can also be designed to help teachers. It is possible to integrate certain repetitive tasks into the games such as the explanation of certain concepts or corrections. This gives the teacher more time to interact with the students, to advise them, and help those in need. The student-teacher relationship is therefore redefined.

In addition, when it comes to training adults, who already have professional experience, games can turn out to be better suited than traditional training (Federation of American Scientist, 2006;[START_REF] Mayo | Games for science and engineering education[END_REF]. This can be explained by the fact that adults find it difficult to "go back to school" and accept criticism from a teacher, who is often younger than them. Learning Games offer an original, less academic way of learning. In addition, the games can be designed so that all forms of judgment come from the game and not the teacher. In Laboratorium of epidemiology9 (Figure 4), the hospital's chief doctor asks the students to write an essay. These essays are actually corrected by the teacher but, the fact that the comments are delivered by the chief doctor gives them more weight and eliminates remonstration [START_REF] Ney | Learning aware environment: a Laboratorium of epidemiological studies[END_REF]. A means to evaluate, track and provide feedback Finally, Learning Games offer the technical means to track the student's actions and automatically evaluate certain skills, a feature which is very important in higher education and professional training [START_REF] Carron | Teaching with Games Based Learning. Management Systems: Exploring and Observing a Pedagogical Dungeon[END_REF]. These tracks can also be analyzed in real time and presented on a monitoring platform for teachers and students themselves. Finally, these tracks can also be used to automatically adapt activities according to the learner's level in order to keep him or her in a state of constant motivation [START_REF] Oostendorp | Adapting the complexity level of a serious game to the proficiency of players[END_REF]. For example, the V3S Virtual Reality training platform offers dynamic situated feedback and progressive learning scenarios by adapting the complexity of the situations to the learner's activity and level [START_REF] Barot | V3S: A Virtual Environment for Risk-Management Training Based on Human-Activity Models[END_REF]. The Reviateck simulators (Figure 5) also analyzes the usage tracks to provide personal feedback on technical and non-technical skills [START_REF] Huguet | Errare Humanum Est: Simulation of communication error among a virtual team in crisis situation[END_REF].

Current Research Questions

Learning Games have much to offer for higher education. However, there remain a number of important issues that need to be addressed before university teachers and professional trainers willingly adopt them. In this section, some of the current limitations and major research questions related to the specificities of Learning Games in higher education will be identified. First of all, research on methods to facilitate the design of custom Learning Games for specific higher education skills will be discussed. Then, the latest research on how to help university teachers resist the skepticism of their students and colleagues when they want to use games in their classes will be presented.

Methods and tools to create custom Learning Games

Specialized websites offer access to Learning Games libraries 10,11,12 . They allow teachers to both find existing Learning Games, they can integrate into their courses, and to share their own Learning Games. These libraries are quite successful with middle school teachers because there is a large body of teachers that are interested in the same learning content, and therefore have the same educational needs. In higher education however, the curriculum is not as precise and teachers create their own pedagogical activities. In this context, it is very difficult to find Learning Games that fit all their requirements.

In order to help teachers create their own custom Learning Games, several researchers have proposed methods and tools such as ScenLRPG [START_REF] Mariais | A Description Grid to Support the Design of Learning Role-Play Games[END_REF], Player-VG [START_REF] Padilla Zea | A Model-Based Approach to Designing Educational Multiplayer Video Games[END_REF] and LEGADEE [START_REF] Marfisi-Schottman | Tools and Methods for Efficiently Designing Serious Games[END_REF] to help them communicate and collaborate with game experts.

Although working with game designers is ideal, it is often too costly and complicated to organize.

Teachers can therefore use Learning Game authoring tools, which allow them to create digital games on their own, without any computer programming skills. Here are a few examples of authoring tools for specific Learning Game types: eAdventure for creating simple point and click games [START_REF] Moreno-Ger | Educational game design for online education[END_REF], Storytec [START_REF] Mehm | User Support in Digital Educational Game Authoring Tools[END_REF] for puzzle games and JEM iNVENTOR [START_REF] Karoui | JEM iNVENTOR: a Mobile Learning Game, Authoring Tool based on a Nested Design Approach[END_REF] for educational mobile treasure hunts. These editors allow teachers to become independent and capable of improving their Learning Games year after year, adapting them to the new curriculum or a modification in the student profile.

An interesting approach, also called game modding, consists in asking players to create new levels for existing games. This movement appeared in the gaming industry with the famous Counter Strike, a "mod" of the game Half Life, designed by two players [START_REF] Kücklich | Precarious Playbour: Modders and the Digital Games Industry[END_REF]. This concept provides several advantages for Learning Games [START_REF] El-Nasr | Learning Through Game Modding[END_REF]. The first advantage is that students are much more likely to create the type of game that they like when following the latest trends (escape games and choose-your-own-adventure type games). The second is the fact that the students will need to acquire a deep understanding of the knowledge and competencies that should be integrated into the Learning Game. Thirdly, teachers then have a large selection of Learning Games to choose from or build on for their future classes. These methods and authoring environments have certainly facilitated the creation of custom Learning Games, but the process remains very complex and time consuming for results that are not always satisfactory. Even though the authoring tools allow teachers (or students) to create their own Learning Games, it is very difficult for them to imagine a gameplay that is fun and that also serves their pedagogical objectives. They might therefore use game mechanics that are not at all adapted to their educational goals. For example, the latest trend is to create Educational Escape Games, even though its core game mechanics (time pressure, mini-puzzles, stress) does not create the right conditions for learning. This type of game, however, is a well-adapted alternative to a quiz in order to verify if the students have understood targeted concepts. The know-how required to select the best game mechanics and create a pedagogically effective Learning Game remains difficult to pinpoint, even for experienced Learning Game designers.

Help teacher convince game skeptics

The negative perception of games in higher education is an important issue. While some university teachers have never been more enthusiastic about Learning Games, the majority remain very skeptical or even openly opposed to such teaching methods, which they see as a waste of time. In this context, it is very difficult for teachers to try using Learning Games, if they are not 100% convinced it is going to work. This implies knowing about game mechanics, the way they can be used to attain educational goals and having several successful examples of Learning Games in mind.

Half a dozen MOOCs on Learning Games and gamification have recently been put online to help teachers acquire such knowledge. For example, MIT offers a MOOC entitled "Design and Development of Games for Learning" 13 , the University of Rotterdam offers a MOOC on Serious Gaming 14 and the University of Pennsylvania offers a MOOC on gamification 15 . Le Mans University 16 , in France, also offers a course to help teachers create their own Learning Games with the help of pedagogical engineers. The student's skepticism toward games can be changed if the teachers clearly explain the educational objectives of the Learning Games and why they have chosen to use a game instead of a classical teaching method. There is a difference between mastering the rules of a game and understanding the educational concepts embedded in it. This a why a phase of reflection and debriefing is necessary after the Learning Game [START_REF] Garris | Games, Motivation, and Learning: A Research and Practice Model[END_REF]. This debriefing serves several purposes such as clarifying the educational purpose of the Learning Game [START_REF] Aldrich | Learning by Doing: A Comprehensive Guide to Simulations, Computer Games, and Pedagogy in e-Learning and Other Educational Experiences[END_REF], fostering 13 https://www.mooc-list.com/course/design-and-development-games-learning-edx 14 https://www.coursera.org/learn/serious-gaming 15 https://www.coursera.org/learn/gamification 16 https://lium.univ-lemans.fr/en/ludifikaction/ metacognition [START_REF] Lederman | Debriefing: Toward a Systematic Assessment of Theory and Practice[END_REF] and facilitating knowledge transfers to other domains and real situations. Some teachers push the concept even further by asking their students to write an "astonishment report" on their experience playing the game and the skills they acquired. This concept is inspired by management methods [START_REF] Vigier | The Astonishment Report: A Pedagogical Tool to Assist Students in Learning from Their International Experience[END_REF]. The evaluation of the students is based on the quality of their essay and not on their actions in the Learning Game. This evaluation method seems effective as it allows the students to learn, without being penalized by errors made in the game, and encourages post analysis and self-refection on the learning outcome.

Educating teachers about game mechanics and how they can be used for teaching and encouraging them to brief/debrief their students are both good ways to help teachers convince Learning Game skeptics among their students and colleagues. However, the best way to convince these skeptics is to provide hard proof that Learning Games enhance learning. The question of evaluation is extremely complicated and has yet to be dealt with. The usual pre and post knowledge tests, carried out with a control group, that does not play the Learning Game, contain too many variables and can easily be flawed.

Cross-references

Augmented Reality and Its Use in Education

Figure 1 .

 1 Figure 1. Train Simulator (Electronic Arts) used for training train drivers (source https://www.origin.com)

Figure 2 .

 2 Figure 2. Starbank The Game (BNP Paribas), used to train bank employees (source http://serious.gameclassification.com)

Figure 3 .

 3 Figure 3. CheckiO used for teaching computer programming (source https://checkio.org/)

Figure 4 .

 4 Figure 4. Laboratorium of Epidemiology (source http://loe.ujf-grenoble.fr)

Figure 5 .

 5 Figure 5. Reviateck simulator for learning technical procedures and social skills (source https://reviatech.com)

Assasin's Creed Ubisoft, https://assassinscreed.ubisoft.com/

https://www.youtube.com/watch?v=TAeut18BHOA

BrainHex questionnaire http://survey.ihobo.com/BrainHex/

The Resuscitation Game (Imaginary), https://www.youtube.com/watch?v=0aDzjJTWUsc&feature=youtu.be

Rail Simulator (Electronic Arts), https://www.origin.com/fra/en-us/store/rail-simulator/rail-simulator

StarBank The Game, http://serious.gameclassification.com/EN/games/14090-Starbank/index.html

Les Aventures de Casey Waren, http://blog.seriousgame.be/les-aventures-de-casey-warren

Laboratorium of Epidemiology, http://loe.ujf-grenoble.fr/content/story


Xavier Cousin

Keywords: Chromatic number, vertices coloring, graph theory Nombre chromatique, coloration de sommets, théorie de graphes

Still other limits of chromatic number

Summary: We present upper bounds for the chromatic number research issued from subgraph of the concerned graph.

Introduction

This paper is to present three limits. We first remind some notions and naming convention. It allows introducing some bounds, maximum limits to the chromatic number based on particular subgraphs. The third limit is the best of the three.

More complex and precise limits can be found in the literature. The goal here is to present my work knowing the third limit can also lower the volume of data manipulated.

The approach progressively starts from a classical limit until an upper bound limiting the research of chromatic number into a subgraph.

Reminders and conventions

• G = (X, A) is the graph which has X as set of vertices and A as set of edges.

• |Y| is the cardinal of a set Y.

• We note d G (x) the degree of x vertex on graph G.

• N is the vertices number of G (or |X|) and M the edges number (or |A|).

• χ is the chromatic number of G corresponding to the minimum number of colours to strongly colour the vertices. In the books of Claude Berge [1] [2], this number is noted  (G).

• A vertex clique is a set of vertices linked together in pairs by edges. We note  the maximum number of a clique.

• We name  the stability number; it is the maximum set of independent vertices (no edge between these vertices two by two). It's the opposite or the complement of the clique.

Graph simplification for the vertices colouring

We will present three maximum limits to the chromatic number χ of any graph G in order to obtain progressively the third which is the optimum. It is based on a subgraph of G.

Theorem 1 -1st bound-(Reminder cf. [1] p326) : χ -1 ≤ max m x  G {d G (x)}. Theorem 2 (Reminder cf. [1] p329) :   x of graph G, χ-critical, d G (x) ≥ χ-1.
Corollary 1:  χ vertices x of any G such that d G (x) ≥ χ -1.

Corollary 2: Let a χ-colouring (S1, ..., Sχ) and di = max m

x  Si {d G (x)} then  i, 1 ≤ i ≤ χ, di ≥ χ-1.

Remark :

• We can sort the stable sets of the χ-colouring by the maximum degree di (by example, d1 ≥ ... ≥ dχ).

• We notice that in the decreasing order case, d1 is the maximum degree in the global graph G.

Xavier Cousin

2/2

Theorem 3 -2nd bound-Let q a maximum number with Y = {x, d G (x) ≥ q-1} and |Y| ≥ q then χ ≤ q.

Proof: (by contradiction) Assume such a set Y with χ = r > q.

Let us take H, χ -critical graph of G, H is a sub-graph of G.

We know by the corollary 2 that  Si, set of the χ-colouring, max m

x  Si (d H (x)) = di ≥ χ -1 then ≥ r-1.

We know also that  r stables, Z1 to Zr then  r vertices x of H such that d H (x) ≥ r-1.

Let W be the set of these vertices, x W and d H (x) ≥ r-1 and |W| ≥ r and r > q.

As H is a sub-graph of G, it is true also for G with d G .

It implies that q is not maximum with such property in G. Thus the contradiction.

Corollary 3 -3rd bound-Let q a maximum number, G1 = (X1, A1), sub-graph of G such that X1 = {x  X, d G1 (x) ≥ q-1} and |X1| ≥ q then χ ≤ q.

The difference with the theorem -2 nd bound-is that the degrees are restricted to the subgraph generated by X1. We can see that the latter is better than the previous one.

Proof: It is identical to the previous one because indeed the degrees of the χ-critical graph are ≥ χ-1 even if they are restricted to G1.

Remark for application in process :

• To find the 2nd bound, it is enough to sort the vertices in decreasing order on the degrees and to constitute vertex by vertex (starting from the highest degree) the set Y wished. Each time, we test if there is a satisfactory value (in the sense of the 2nd bound) for q. The first value found concerns the first highest degrees and is therefore maximum (the other values will be based on lower degrees).

• For the 3 rd bound, the process is as follows. We loop of q = | X | (or degree max m + 1) down to q = 1 by constructing a set Y (d G (x) ≥ q-1). Y is disconnected from the graph G. The other vertices of G are not affected by this value of q. We recursively purify Y by removing the vertices whose degree becomes < q-1. The process stops as soon as there is nothing left to clean. This means that all the vertices of the purified set are of degree (in the set Y) ≥ q-1. There are two possibilities. Let the number of vertices ≥ q-1 be strictly less than q then we decrement the value of q and start again. In the other case, |Y| ≥ q. The other values of q having failed, and the process being decreasing, this value is therefore maximal.

Corollary 4 :   vertices x of any G such that d G (x) ≥  -1.

Corollary 5 :   vertices x of any G = (X, A) such that d G (x) < |X| - + 1.

Indeed, by taking the complementary graph , to search

 corresponds to search  in G. Then     vertices x of such that d G (x) ≥    -1. But d G (x)= |X| -d G (x). Then   G vertices x of G such that: |X| -d G (x) ≥  G -1 or by multiplying by -1 : d G (x) -|X| < 1 - G or d G (x) < |X| - G + 1.

Conclusion and synthesis

In conclusion, to adopt the process of the third bound is good enough to obtain the chromatic number by reducing the graph to be coloured.

In summary we have a bound greater than χ: let q maximum number, SG = (SX, SA), sub-graph of G such that SX = {x X, d SG (x) ≥ q-1} and |SX| ≥ q then χ ≤ q.

[1] C. Berge, Graphes, Bordas, Paris, 1983.

[2] C. Berge, Hypergraphes, Combinatoire des ensembles finis, Bordas, Paris, 1987.


domain: spi phys math info

text:
Ye Shi
email: ye.shi@student.uts.edu.au

Hoang Duong
email: tuan.hoang@uts.edu.au

Pierre Apkarian
email: pierre.apkarian@onera.fr

Andrey V Savkin
email: a.savkin@unsw.edu.au

H D Tuan

Global Optimal Power Flow over Large-Scale Power Transmission Networks

Keywords: Optimal power flow (OPF) problem, large-scale transmission networks, rank-one matrix constraint, nonsmooth optimization, semi-definite programming (SDP)

  L'archive ouverte pluridisciplinaire

Introduction

Smart grids are operated by the advanced distribution management system (DMS), which is responsible for supervisory control and data acquisition in reactive dispatch, voltage regulation, contingency analysis, capability maximization and other smart operations. The optimal power flow (OPF) problem, which determines a steady state operating point that minimizes the cost of electric power generation or the transmission loss is the backbone of DMS (see e.g. [START_REF] Carpentier | Contribution to the economic dispatch problem[END_REF][START_REF] Huneault | A survey of the optimal power flow literature[END_REF][START_REF] Momoh | A review of selected optimal power flow literature to 1993. part I: Nonlinear and quadratic programming approaches[END_REF][START_REF] Pandya | A survey of optimal power flow methods[END_REF] and references therein). Mathematically, the OPF problem is highly nonlinear and nonconvex due to numerous quadratic equality and indefinite quadratic inequality constraints for bus interconnections, hardware operating capacity and the balance between power demand and supply. These nonlinear constraints are mathematically troublesome so the state-ofthe-art nonlinear optimization solvers may converge to just stationary points (see [START_REF] Bukhsh | Local solutions of the optimal power flow problem[END_REF] and references therein), which are not necessarily feasible. To handle these nonlinear constraints, it is common to reformulate them as linear constraints on the outer product W = V V H of the voltage vector V = (V 1 , V 2 , ..., V n ) T ∈ C n . As a result, the OPF problem is recast by a semi-definite program (SDP) plus the additional rank-one constraint on outer product matrix W [START_REF] Bai | Semidefinite programming for optimal power flow problems[END_REF][START_REF] Lavaei | Zero duality gap in optimal power flow problem[END_REF]. The rank-one constraint is then dropped for semi-definite relaxation (SDR). However, the optimal solution of SDR is of rank-more-than-one in general and cannot help retrieval of a feasible point or stationary point of the OPF problem [START_REF] Lavaei | Zero duality gap in optimal power flow problem[END_REF][START_REF] Madani | Convex relaxation for optimal power flow problem: Mesh networks[END_REF][START_REF] Madani | Promises of conic relaxation for contigency-constrained optimal power flow problem[END_REF]. In [START_REF] Shi | Nonsmooth optimization for optimal power flow over transmission networks[END_REF], we have extended the technique of [START_REF] Phan | Nonsmooth optimization for efficient beamforming in cognitive radio multicast transmission[END_REF] for solving the beamforming optimization problems in signal processing to optimize the outer product matrix W , which works very well and is practical for moderate-scale power distribution networks up to n = 300 buses. There is another approach (see e.g. [START_REF] Josz | Moment/sum-of-squares hierarchy for complex polynomial optimization[END_REF] and references therein), which is based on hierarchies of moment-based relaxation for nonconvex quadratic problems to tackle large networks with simple nonconvex constraints. Power transmission networks in modern smart grids are often devised with a few thousand buses [START_REF] Bamberger | Vision and Strategy for Europes Electricity Networks of the Future[END_REF][START_REF] Hallberg | Smart Grids and Networks of the Future -EURELECTRIC Views[END_REF][START_REF] Hashmi | Survey of smart grid concepts, architectures, and technological demonstrations worldwide[END_REF]. Under a such large number n of buses it is impossible to use the single matrix W ∈ C n×n , which involves n(n+1)/2 ≈ O(10 7 ) complex variables. On the other hand, the number of the flow lines for bus connection is relatively moderate so only a small portion of the crossed nonlinear terms V k V * m appears in the nonlinear constraints. The common approach is to use the outer products of overlapped groups of the voltage variables to cover them [START_REF] Molzahn | Implementation of a large-scale optimal power flow solver based on semidefinite programming[END_REF][START_REF] Andersen | Reducedcomplexity semidefinite relaxations of optimal power flow problems[END_REF][START_REF] Madani | Promises of conic relaxation for contigency-constrained optimal power flow problem[END_REF]. All rank-one constraints on these outer products are then dropped for SDR. Obviously, the optimal solution of this SDR usually is not of rank-one and thus does not have any physical meaning. There is no technique to retrieve a feasi-ble rank-one point from the rank-more-than-one solution of SDR. 1 Multiple matrix rank constrained optimization has received a great attention due to its potential application in robust control synthesis [START_REF] Apkarian | A sequential SDP/Gauss-Newton algorithm for rank-constrained LMI problems[END_REF][START_REF] Apkarian | A trust region spectral bundle method for nonconvex eigenvalue optimization[END_REF] but to our best knowledge there is no effective computation so far. The contribution of this paper is two-fold:

• An effective decomposition for large-scale OPF problems, which involves essentially reduced numbers of the rank-one constraints on matrices of moderate size for expressing the network nonlinear constraints;

• A new iterative procedure for rank-one constrained optimization, which is practical for computational solutions of large-scale indefinite quadratic programming. Simulations for the large-scale OPF problems show that it is capable of finding the global optimal solution with the computational tolerance less than 0.01%.

The paper is structured as follows. Section 2 is devoted to the OPF problem formulation and its difficulties. Its computational solution is developed in Section 3. Section 4 provides simulation to show the efficiency of our method.

The conclusions are drawn in Section 5. Notation. j denotes the imaginary unit; M 0 means that M is a Hermitian symmetric positive semi-definite matrix; rank(M ) is the rank of the matrix M ; (•) and (•) denote the real and imaginary parts of a complex quantity; a ≤ b for two complex numbers a and b is componentwise understood, i.e. (a) ≤ (b) and (a) ≤ (b); ., . is the dot product of matrices, while diag{A i } denotes the matrix with diagonal blocks A i and zero off-diagonal blocks; the cardinality of a set L is denoted by |L|.

Optimal power flow problem and challenges

Consider an AC electricity transmission network with a set of n buses N := {1, 2, • • • , n}. The buses are connected through a set of flow lines L ⊆ N × N , i.e. bus m is connected to bus k if and only if (m, k) ∈ L. Accordingly, N (k) := {m ∈ N : (m, k) ∈ L}. The power demanded at bus k ∈ N is S L k = P L k + jQ L k , where P L k and Q L k are the real and reactive power. A subset G ⊆ N of buses is supposed to be connected to generators. Any bus k ∈ N \G is thus not connected to generators.

Other physical parameters are following [START_REF] Hadi | Power system analysis[END_REF][START_REF] Dommel | Optimal power flow solutions[END_REF][START_REF] Zimmerman | MAT-POWER: steady-state operations, planning, and analysis tools for power systems research and education[END_REF]:

• Y = [y km ] (k,m)∈N ×N ∈ C n×n
is the admittance matrix [START_REF] Zimmerman | MAT-POWER: steady-state operations, planning, and analysis tools for power systems research and education[END_REF]. Each y km is the mutual admittance between bus k and bus m, so y km = y mk ∀ (k, m) ∈ L.

• V is the complex voltage vector,

V = [V 1 , V 2 , • • • , V n ] T ∈ C n
, where V k is the complex voltage injected to bus k ∈ N .

1 There is an algorithm of finding a rank-one solution [18, Alg. 1], which however is applicable to simple nonconvex constraints and is not guaranteed to convergence

• I is the complex current vector, I = Y V = [I 1 , I 2 , • • • I n ] T ∈ C n
, where I k is the complex current injected to bus k ∈ N .

• I km is the complex current in the power line (k, m) ∈ L, m∈N (k)

I km = I k = m∈N (k)
y km V m .

• S km = P km +jQ km is the complex power transferred from bus k to bus m, where P km and Q km represent the real and reactive transferred power.

• S G k = P G k + jQ G k is the complex power injected by bus k ∈ G, where P G k and Q G k represent the real and reactive generated power.

For each bus k, it is obvious that

S G k -S L k = (P G k -P L k ) + j(Q G k -Q L k ) = V k I * k = V k m∈N (k) V * m y * km .
Therefore, the real generated power P G k and reactive generated power Q G k at bus k are the following nonconvex quadratic functions of the bus voltage vector variable

V := (V 1 , V 2 , ..., V n ) T ∈ C n : P G k = P L k + ( m∈N (k) V k V * m y * km )
and For illustrative purpose, Figure 1 provides a diagram of WB5-Five Bus Network [START_REF] Bukhsh | Test case archive of optimal power flow (OPF) problems with local optima[END_REF] with N = {1, 2, 3, 4, 5}. It is seen from the bus connection in this Figure that N (1) = {2, 3} because buses 2 and 3 are connected to bus 1, N (2) = {1, 3, 4}, N (3) = {1, 5}, N (4) = {2, 5} and N (5) = {3, 4}. Also G = {1, 5} because buses 1 and 5 are connected to generators.

Q G k = Q L k + ( m∈N (k) V k V * m y * km ).
The objective of OPF is to minimize either the total cost of the generated power defined by [START_REF] Dommel | Optimal power flow solutions[END_REF]

f (P G ) = k∈G (c k2 P 2 G k + c k1 P G k + c k0 ), (1)
with given c k2 > 0, c k1 and c k0 and real active generated power P G , or the total transmission losses defined by

F loss (P G ) = k∈G (P G k -P L k ) - k∈N \G P L k . (2)
In this paper, we consider only the objective function [START_REF] Carpentier | Contribution to the economic dispatch problem[END_REF].

It is quite clear that our approach can be adjusted to solve the objective (2) in a straightforward manner. Write the objective (1) as the following function of the bus voltages V :

f (V ) = k∈G [c k2 (P L k + ( m∈N (k) V k V * m y * km )) 2 +c k1 (P L k + ( m∈N (k) V k V * m y * km )) + c k0 ]. (3)
Accordingly, the OPF problem is formulated as follows

min V ∈C n f (V ) s.t. (4a) -P L k -jQ L k = m∈N (k) V k V * m y * km , k ∈ N \ G, (4b)
P min G k ≤ P L k + ( m∈N (k) V k V * m y * km ) ≤ P max G k , k ∈ G (4c) Q min G k ≤ Q L k + ( m∈N (k) V k V * m y * km ) ≤ Q max G k , k ∈ G (4d) V min k ≤ |V k | ≤ V max k , k ∈ N , ( 4e
)
|S km | = |V k V * m y * km | ≤ S max km , ∀(k, m) ∈ L (4f) |V k -V m | ≤ V max km , (k, m) ∈ L, (4g) |arg(V k ) -arg(V m )| ≤ θ max km , (k, m) ∈ L, (4h)
where (4b) is the equation of the balance between the demand and supply power at bus k ∈ N \ G, (4c)-(4d) are the power generation bounds with the lower bound

P min G k , Q min G k and upper bound P max G k , Q max G k of
the real power reactive power generations, (4e) are the voltage amplitude bounds, (4f)-(4h) are capacity limitations with the line currents between the connected buses constrained by (4f), while (4g)-(4h) guarantee the voltage balance in terms of their magnitude and phases [START_REF] Zimmerman | MAT-POWER: steady-state operations, planning, and analysis tools for power systems research and education[END_REF]. One can see that (4) is a highly nonlinear optimization subject to nonlinear inequality and equality constraints and thus poses a real challenge for its computation. A common approach is to use the slack variables

W km = V k V * m , k = 1, ..., n; m = 1, ..., n (5)
and then recast the problem (4) in V ∈ C n to the following

problem in W = [W km ] k,m=1,...,n ∈ C n×n : min W ∈C n×n F (W ) s.t. W 0, (6a) -P L k -jQ L k = m∈N (k) W km y * km , k ∈ N \ G, (6b)
P min G k ≤ P L k + ( m∈N (k) W km y * km ) ≤ P max G k , (6c)
k ∈ G, Q min G k ≤ Q L k + ( m∈N (k) W km y * km ) ≤ Q max G k , (6d)
k ∈ G, (V min k ) 2 ≤ W kk ≤ (V max k ) 2 , k ∈ N , (6e)
|W km y * km | ≤ S max km , (k, m) ∈ L, (6f)
W kk + W mm -W km -W mk ≤ (V max km ) 2 , ( 6g
) (k, m) ∈ L, (W km ) ≤ (W km ) tan θ max km , (k, m) ∈ L, (6h) rank(W ) = 1, (6i)
where

F (W ) = k∈G [c k2 (P L k + ( m∈N (k) W km y * km )) 2 + c k1 (P L k + ( m∈N (k) W km y * km )) + c k0 ],
which is convex quadratic in W km , while all constraints (6b)-(6h) are linear. All the problem nonconvexity is now concentrated in the single rank-one constraint (8b), which is then dropped for SDR. If the solution of this SDR is of rank-one, then it obviously leads to the global solution of the nonconvex optimization problem [START_REF] Bai | Semidefinite programming for optimal power flow problems[END_REF]. Otherwise even a feasible point of ( 6) is hardly retrieved from the SDR solution. Moreover formulation (6) cannot be practically used for n more than a few thousands in large-scale networks as it involves n(n + 1)/2 variables-a prohibitively large number. The next section will address and resolve all these issues.

Nonsmooth optimization based solution

The first issue is to decompose large-size matrix W in (6) into matrices of smaller size to make the problem tractable even with limited computational power. This is also prompted by the fact that there is only a small portion of the crossed terms V k V * m appearing in the nonlinear constraints (4b)-(4h) so the large-size matrix variable W ∈ C n×n contains many redundant terms V k V * m . The main result of [START_REF] Molzahn | Implementation of a large-scale optimal power flow solver based on semidefinite programming[END_REF][START_REF] Andersen | Reducedcomplexity semidefinite relaxations of optimal power flow problems[END_REF][START_REF] Madani | Promises of conic relaxation for contigency-constrained optimal power flow problem[END_REF] is to decompose the set N := {1, 2, ..., n} of buses into I overlapped subsets N i = {i 1 , ..., i Ni } of buses, called bags, such that i ∈ N (i +1 ), = 1, ..., N i -1 and i Ni ∈ N (i 1 ), for each i = 1, 2, ..., I, i.e. the buses in the same bag are serially connected. The set of bags can be reset to make bags of relatively same size. Define the Hermitian symmetric matrix variables [START_REF] Bai | Semidefinite programming for optimal power flow problems[END_REF] we have the following equivalent reformulation for [START_REF] Bai | Semidefinite programming for optimal power flow problems[END_REF] min

W i = [W i k im ] k,m=1,..,Ni ∈ C Ni×Ni , i = 1, 2, ..., I. (7) By replacing W km = V k V * m in
W =diag{W i } F (W ) s.t. (6b) -(6h), (8a)
W i 0, i = 1, ..., I, (8b) rank(W i ) = 1, i = 1, ..., I. (8c)
Reference [START_REF] Molzahn | Implementation of a large-scale optimal power flow solver based on semidefinite programming[END_REF][START_REF] Andersen | Reducedcomplexity semidefinite relaxations of optimal power flow problems[END_REF] just dropped all rank-one constraints in (8c) for SDR without any justification. Reference [START_REF] Madani | Promises of conic relaxation for contigency-constrained optimal power flow problem[END_REF] also dropped all rank-constraints in (8c) but then used a penalized SDR for locating low-rank semi-definite matrices W i in [START_REF] Lavaei | Zero duality gap in optimal power flow problem[END_REF]. Based on these low-rank matrices, [9, Sec. IV] also proposed to find rank-one matrices, which however are not necessarily feasible to [START_REF] Madani | Convex relaxation for optimal power flow problem: Mesh networks[END_REF].

The variable number in [START_REF] Madani | Convex relaxation for optimal power flow problem: Mesh networks[END_REF] is

I i=1 N i (N i + 1)/2.
To keep this number reasonably moderate, it is desired that both I and N i are sufficiently moderate. However, one can see that the above described decomposition [START_REF] Molzahn | Implementation of a large-scale optimal power flow solver based on semidefinite programming[END_REF][START_REF] Andersen | Reducedcomplexity semidefinite relaxations of optimal power flow problems[END_REF][START_REF] Madani | Promises of conic relaxation for contigency-constrained optimal power flow problem[END_REF] leads to a large number I of bags as well as few large size N i that result in many rank-one constraints in (8c), which are much less probably satisfied by solving SDR.

Our first step toward to computation of ( 4) is to develop a new decomposition with many fewer bags involved. Recalling that N (k) is the set of the buses that are connected to bus k, the cardinality |N (k)| is small in largescale networks. We resort N = {1, 2..., n} as

N = {N 1 , ..., N n } such that the cardinality |N (N k )| is in decreased or- der: |N (N 1 )| ≥ |N (N 2 )| ≥ ... ≥ |N (N n )|.
Accordingly, the first bag of buses is defined as N 1 = N (N 1 ). The second bag is defined as

N 2 = {i ∈ N (N 2 ) : {i, N 2 } ⊂ N 1 }.
Note that the crossed term V i V * N2 is already treated in the previous bag N 1 whenever {i, N 2 } ⊂ N 1 so we exclude such bus i in defining bag N 2 . Similarly, for ≥ 3 the -th bag is defined as

N = {i ∈ N (N ) : {i, N } ⊂ N ∀1 ≤ ≤ -1}
to exclude those buses i, whose crossed term V i V * N already is treated in a previous bag. As each |N i | is obviously small, such decomposition is very efficient, leading to a substantial reduction of involved bags in comparison with that used in [START_REF] Molzahn | Implementation of a large-scale optimal power flow solver based on semidefinite programming[END_REF][START_REF] Andersen | Reducedcomplexity semidefinite relaxations of optimal power flow problems[END_REF][START_REF] Madani | Promises of conic relaxation for contigency-constrained optimal power flow problem[END_REF]. Our next step is to tackle the numerous difficult rankone constraints in (8c), not dropping them for SDR as in all the previous works. Firstly we express I rank-one constraints in (8c) by the following single spectral constraint

I i=1 (Trace(W i ) -λ max (W i )) = 0, (9)
where λ max (W i ) stands for the maximal eigenvalue of W (i) .

Indeed, (8b) implies Trace(W i ) -λ max (W i ) ≥ 0 ∀ i, so (9) means that Trace(W i ) = λ max (W i ), i.e. W i has only one nonzero eigenvalue so it is of rank-one. The nonnegative quantity I i=1 (Trace(W i ) -λ max (W i )) can therefore be used to measure the degree of satisfaction of the rankone constraints [START_REF] Madani | Promises of conic relaxation for contigency-constrained optimal power flow problem[END_REF]. Without squaring, the penalization I i=1 (Trace(W i ) -λ max (W i )) is exact, meaning that (9) can be satisfied by a minimizer of the problem min 8b), [START_REF] Shi | Nonsmooth optimization for optimal power flow over transmission networks[END_REF] with a finite value of µ > 0 (see e.g. [START_REF] Bonnans | Numerical Optimization Theoretical and Practical Aspects[END_REF]Chapter 16]). This is generally considered as a sufficiently nice property to make such exact penalization attractive. For any W i,(κ) feasible for the convex constraints (6b)-(6h), (8b), function λ max (W i ) is nonsmooth and is lower bounded by

W =diag{W i } F µ (W ) := F (W ) + µ I i=1 (Trace(W i ) -λ max (W i )) s.t. (6b) -(6h), (
λ max (W i ) = max ||w||=1 w H W i w ≥ (w i,(κ) max ) H W i w i,(κ) max , (11)
where w i,(κ) max is the normalized eigenvector corresponding to the eigenvalue λ max (W i,(κ) ), i.e.

λ max (W i,(κ) ) = (w i,(κ) max ) H W i,(κ) w i,(κ) max . (12)
Accordingly, µλ max

(W i )-µλ max (W i,(κ) ) ≥ µ((w i,(κ) max ) H W i .w i,(κ) max -(w i,(κ) max ) H W i,(κ) w i,(κ) max ) = µ w i,(κ) max (w i,(κ) max ) H , W i - W i,(κ) , so µw i,(κ) max (w i,(κ) max ) H is a subgradient of the function µλ max (W i ) at W i,(κ) . Then µdiag{w i,(κ) max (w i,(κ) max ) H } is a subgradient of the function µ I i=1 λ max (W i ) at diag{W i,(κ) }.
The following SDP provides an upper bound for the nonconvex optimization problem [START_REF] Shi | Nonsmooth optimization for optimal power flow over transmission networks[END_REF] min [START_REF] Phan | Nonsmooth optimization for efficient beamforming in cognitive radio multicast transmission[END_REF]. Suppose that W (κ+1) = diag{W i,(κ+1) } is the optimal solution of SDP [START_REF] Bamberger | Vision and Strategy for Europes Electricity Networks of the Future[END_REF]. Since W (κ) = diag{W i,(κ) } is also feasible to [START_REF] Bamberger | Vision and Strategy for Europes Electricity Networks of the Future[END_REF] with

W =diag{W i } F (κ) (W ) := F (W ) + µ I i=1 (Trace(W i ) -(w i,(κ) max ) H W i w i,(κ) max ) s.t. (6b) -(6h), (8b) (13) because F (κ) (diag{W i }) ≥ F µ (diag{W i }) ∀ W i 0 ac- cording to
F µ (W (κ) ) = F (κ) (W (κ) ), it is true that F µ (W (κ+1) ) ≤ F (κ) (W (κ+1) ) ≤ F (κ) (W (κ) ) = F µ (W (κ) ), so W (κ+1
) is a better feasible point of (10) than W (κ) . Initialized by any feasible point W (0) = diag{W i,(0) } of SDP constraint (8a)-(8b), the sequence {W (κ) } = {diag{W i,(κ) }} with W (κ+1) = diag{W i,(κ+1) } iteratively generated as the optimal solution of SDP ( 13) is a sequence of improved feasible points of the nonconvex optimization problem [START_REF] Shi | Nonsmooth optimization for optimal power flow over transmission networks[END_REF].

Since W (κ) are uniformly bounded, the sequence {W (κ) } has a limit point W = diag{ W i }, which is the optimal solution of the optimization problem min

W =diag{W i } F (W ) + µ I i=1 (Trace(W i ) -( wi max ) H W i wi max ) s.t. (6b) -(6h), (8b), ( 14
)
where wi max is the normalized eigenvector corresponding to the eigenvalue λ max ( W i ) of W i . Particularly,

F (W ) + µ I i=1 (Trace(W i ) -( wi max ) H W i wi max ) ≥ F (W ) + µ I i=1 (Trace( W i ) -( wi max ) H W i wi max ),
or equivalently, under the definition g(W

) = F (W ) + µ I i=1 Trace(W i ), g(W ) -g( W ) -µdiag{ wi max ( wi max ) H }, W -W ≥ 0 for all feasible points W = diag{W i } in (6b)-(6h), (8b).
As a result, W is the optimal solution of the convex optimization problem min 6b) -(6h), (8b), so it must satisfy the optimality condition ∇g( W ) -µdiag{ wi max ( wi max ) H }, W -W ≥ 0 for all feasible points diag{W i } in (6b)-(6h), (8b). The latter is also the first order necessary optimality condition for [START_REF] Shi | Nonsmooth optimization for optimal power flow over transmission networks[END_REF] because µdiag{ wi max ( wi max ) H } is a subgradient of the function µ I i=1 λ max (W i ) at W . As our simulations will show, W is indeed the global optimal solution of ( 10) and [START_REF] Madani | Convex relaxation for optimal power flow problem: Mesh networks[END_REF].

W =diag{W i } g(W ) -µdiag{ wi max ( wi max ) H }, W -W s.t. (
However, unlike [START_REF] Shi | Nonsmooth optimization for optimal power flow over transmission networks[END_REF] with only a single rank-one constrained matrix, although quantity

I i=1 (Trace(W i,(κ) ) -λ max (W i,(κ) )) ( 15
)
in ( 10) is iteratively decreased, not all individual quantities

Trace(W i,(κ) ) -λ max (W i,(κ) ) ( 16)

are iteratively decreased so the rank of each matrix W i,(κ) is no longer iteratively reduced to one as expected. Worse, W i,(κ) is rank-one but the rank of W i,(κ+1) in the next iteration may turn to be more than one with

Trace(W i,(κ+1) ) -λ max (W i,(κ+1) ) > Trace(W i,(κ) ) -λ max (W i,(κ) ).
Consequently, it is very difficult to achieve rank-one for all W i,(κ) as desired. It is also impossible to add a "weight" to each term under the sum in the objective in [START_REF] Bamberger | Vision and Strategy for Europes Electricity Networks of the Future[END_REF] to handle the individual convergence of Trace(W i ) -λ max (W i ).

We now develop a systematic way to resolve this issue as follows. For κ = 0, 1, ..., and W (κ) = diag{W (i,(κ)) } define

L (κ) = {i ∈ {1, ..., I} : rank(W i,(κ) ) = 1} (17)
and generate W (κ+1) = diag{W i,(κ+1) } as the optimal solution of the following SDP instead of SDP [START_REF] Bamberger | Vision and Strategy for Europes Electricity Networks of the Future[END_REF] min

W =diag{W i } F (W ) + µ I i=1 [Trace(W i ) -(w i,(κ) max ) H W i w i,(κ) max ] s.t. (6b) -(6h), (8b), (18a) Trace(W i ) -(w i,(κ) max ) H W i w i,(κ) max ≤ tol , i ∈ L (κ) .( 18b
)
Note that Trace(W i ) ≥ w H W w for all ||w|| = 1 and it is obvious that rank(W i ) = 1 if and only if Trace(W i )w H max W i w max = 0 for some normalized w max . Therefore, the constraint (18b) for some tolerance tol is introduced to warrant the rank-one of all W i,(κ+1) , i ∈ L (κ) . As a result L (κ) ⊂ L (κ+1) and L (κ) → {1, ..., I} is expected to have all W i,(κ) of rank-one. Unlike (13), the iterations (18) leads to achieving rank-one of all W i while the objective function F µ is still decreased.

In summary, we propose the following Large-Scale Nonsmooth Optimization Algorithm (Large-scale NOA) for the multiple rank-one constrained optimization problem [START_REF] Shi | Nonsmooth optimization for optimal power flow over transmission networks[END_REF].

Initialization. Solve SDP min W =diag{W i } F (W ) s.t. (6b) -(6h), (8b) (19)
to generate W (0) := diag{W i,(0) }. If rank(W i,(0) ) ≡ 1 stop: W (0) is the global solution of the nonconvex optimization problem [START_REF] Madani | Convex relaxation for optimal power flow problem: Mesh networks[END_REF]. Otherwise set κ = 0 and define L (κ) by ( 17). κ-th iteration. For κ = 0, 1, .., solve [START_REF] Molzahn | A Laplacian-based approach for finding near globally optimal solutions to OPF problems[END_REF] to generate W (κ+1) := diag{W i,(κ+1) }. Reset κ = κ + 1 and define L (κ) by [START_REF] Andersen | Reducedcomplexity semidefinite relaxations of optimal power flow problems[END_REF]. Stop whenever L (κ) = {1, ..., I}. Otherwise go to the next iteration.

Simulation results

The computation facilities for our implementation are Processor Intel(R) Core i5-3470 CPU @3.20GHz, Matlab version R2013b and CVX with SDPT3. We set the tolerances = tol = 10 -5 and the penalty parameter µ = 10 6 which makes the penalty term µ I i=1 (Trace(W (i,(0)) )λ max (W i,(0) )) at similar magnitude with the objective F (W (0) ). The data source for all examples is Matpower version 5.1 [START_REF] Zimmerman | MAT-POWER: steady-state operations, planning, and analysis tools for power systems research and education[END_REF]. All examples were considered in [START_REF] Molzahn | Implementation of a large-scale optimal power flow solver based on semidefinite programming[END_REF][START_REF] Andersen | Reducedcomplexity semidefinite relaxations of optimal power flow problems[END_REF][START_REF] Madani | Promises of conic relaxation for contigency-constrained optimal power flow problem[END_REF] by SDR only. We recall that T is the number of matrix variables W i in the OPF problem [START_REF] Madani | Convex relaxation for optimal power flow problem: Mesh networks[END_REF] and L (κ) is defined by [START_REF] Andersen | Reducedcomplexity semidefinite relaxations of optimal power flow problems[END_REF] is the set of rank-one matrices W i found after κ-th iteration. The capability of our large-scale NOA in locating the global optimal solution of the OPF problems is demonstrated by showing that the global optimality tolerance (GOT) of its found solution defined as the found value-lower bound lower bound is almost zero.

The numerical examples are presented as follows.

Polish-2383wp system

There are n = 2383 buses, 327 generators and 2896 transmission lines, leading to 2056 nonlinear constraints in (4b). Initialization. A lower bound 1.8490×10 6 of ( 8) is found by solving SDP [START_REF] Apkarian | A sequential SDP/Gauss-Newton algorithm for rank-constrained LMI problems[END_REF]. |L (0) | = 1210 and there are 32 matrices W i,(0) of rank-more-than-one. Their largest size (smallest size, resp.) is 10 × 10 (2 × 2, resp.). Stage 1. |L (10) | = 1234 is achieved. There are 8 matrices W i, (10) of rank-more-than-one. Their largest size (smallest size, resp.) is 10 × 10 and (3 × 3, resp.). Stage 2. |L (19) | = 1237 is achieved. There are 5 matrices W i, (20) of rank-more-than-one. Their largest size (smallest size, resp.) is 9 × 9 (3 × 3, resp.). Stage 3. |L (25) | = I = 1242 is achieved. The found value of the objective is 1.8408 × 10 6 with GOT 4.3267e -04.

Polish-2736sp system

There are n = 2736 buses, 420 generators and 3504 transmission lines, which lead to 2316 nonlinear constraints in (4b). Initialization. A lower bound 1.3041 × 10 6 of ( 8) is obtained by solving SDP [START_REF] Apkarian | A sequential SDP/Gauss-Newton algorithm for rank-constrained LMI problems[END_REF]. |L (0) | = 1534 and there are 4 matrices W i,(0) of rank-more-than-one. Their largest size (smallest size, resp.) is 6 × 6 (4 × 4, resp.). Stage 1. |L (9) | = I = 1538 is achieved. The found value of the objective is 1.3042 × 10 6 with GOT 7.6681e -05.

Polish-2737sop system

There are n = 2737 buses, 399 generators and 3506 transmission lines, which lead to 2338 nonlinear constraints in (4b). Initialization. A lower bound 7.7571 × 10 5 of ( 8) is obtained by solving SDP [START_REF] Apkarian | A sequential SDP/Gauss-Newton algorithm for rank-constrained LMI problems[END_REF]. |L (0) | = 1532 and there are 6 matrices W i,(0) of rank-more-than-one. Their largest size (smallest size, resp.) is 6 × 6 (3 × 3, resp.). Stage 1. |L (2) | = I = 1538 is achieved. The found value of the objective is 7.7572 × 10 5 with GOT 1.2891e -05.

Polish-2746wop system

There are n = 2746 buses, 514 generators and 3514 transmission lines, which lead to 2232 nonlinear constraints in (4b). Initialization. A lower bound 1.2039 × 10 6 of ( 8) is obtained by solving SDP [START_REF] Apkarian | A sequential SDP/Gauss-Newton algorithm for rank-constrained LMI problems[END_REF]. |L (0) | = 1538 and there are 8 matrices W i,(0) of rank-more-than-one. Their largest size (smallest size, resp.) is 6 × 6 (3 × 3, resp.). Stage 1. |L (2) | = I = 1546 is achieved. The found value of the objective is 1.2040 × 10 6 with GOT 8.3063e -05.

Polish-2746wp system

There are n = 2746 buses, 520 generators and 3514 transmission lines, which lead to 2226 nonlinear constraints in (4b). Initialization. A lower bound 1.626590 × 10 6 of ( 8) is obtained by solving SDP [START_REF] Apkarian | A sequential SDP/Gauss-Newton algorithm for rank-constrained LMI problems[END_REF]. |L (0) | = 1545 and there are 2 matrices W i,(0) of rank-more-than-one. Their size is 4 × 4. Stage 1. |L (1) | = I = 1547 is achieved. The found value of the objective is 1.626591 × 10 6 with GOT 6.1478e -07.

Polish-3012wp system

There are n = 3012 buses, 502 generators and 3572 transmission lines, which lead to 2510 nonlinear constraints in (4b). Initialization. A lower bound 2.5717 × 10 6 of ( 8) is obtained by solving SDP [START_REF] Apkarian | A sequential SDP/Gauss-Newton algorithm for rank-constrained LMI problems[END_REF]. |L (0) | = 1682 and there are 7 matrices W i,(0) of rank-more-than-one. Their largest size (smallest size, resp.) is 7 × 7 (2 × 2, resp.). Stage 1. |L (4) | = I = 1689 is achieved. The found value of the objective is 2.5727 × 10 6 with GOT 3.8885e -04.

Polish-3120sp system

There are n = 3120 buses, 505 generators and 3693 transmission lines, which lead to 2615 nonlinear constraints in (4b). Initialization. A lower bound 2.1314 × 10 6 of ( 8) is obtained by solving SDP [START_REF] Apkarian | A sequential SDP/Gauss-Newton algorithm for rank-constrained LMI problems[END_REF]. |L (0) | = 1749 and there are 8 matrices W i,(0) of rank-more-than-one. Their largest size (smallest size, resp.) is 8 × 8 (2 × 2, resp.). Stage 1. |L (9) | = I = 1757 is achieved. The found value of the objective is 2.1391 × 10 6 with GOT 0.0036.

Numerical summary

One can observe that GOT of the solutions computed by the large-scale NOA is very small, proving its capability to provide the global solution of [START_REF] Madani | Convex relaxation for optimal power flow problem: Mesh networks[END_REF]. Table 1 and Table 2 summarize the main points in our simulation. The second and third columns of Table 1 are the number I of bags in (7) by our decomposition and by that in [START_REF] Madani | Promises of conic relaxation for contigency-constrained optimal power flow problem[END_REF], while the fourth and fifth columns give the maximum size N i in [START_REF] Lavaei | Zero duality gap in optimal power flow problem[END_REF]. One can see that both I and the maximum N i by our decomposition are substantially smaller than their counterparts by [START_REF] Madani | Promises of conic relaxation for contigency-constrained optimal power flow problem[END_REF]. This leads to far smaller numbers of variables in [START_REF] Madani | Convex relaxation for optimal power flow problem: Mesh networks[END_REF], which are provided in the seventh and eighth columns. The number n(n + 1)/2 of complex variables in ( 6) is also provided in the sixth column to contrast to the number of complex variables in [START_REF] Madani | Convex relaxation for optimal power flow problem: Mesh networks[END_REF] in the seventh column. Furthermore, the second column of Table 2 provides the best values of (4) found by our large-scale NOA, which are far smaller than ones in the third and fourth columns found by [START_REF] Madani | Promises of conic relaxation for contigency-constrained optimal power flow problem[END_REF] and Matpower6.0 [START_REF] Zimmerman | MAT-POWER: steady-state operations, planning, and analysis tools for power systems research and education[END_REF] (using an interior point method), respectively. In short, our computation approach to the OPF problem (8) outperforms other existing approaches in terms of computational efficiency and performance.

Conclusion

The OPF problems over power transmission network are large-scale optimization problems, which involve a large number of quadratic equality and indefinite quadratic inequality constraints and thus are difficult computationally. We have developed a large-scale nonsmooth optimization algorithm to compute their optimal solutions, which is efficient and practical for large-scale power transmission networks of a few thousands of buses. Applications of the developed large-scale NOA to the OPF problems over three-phase power transmission networks are currently under investigation.

Figure 1 :

 1 Figure1: WB5-Five Bus Network[START_REF] Bukhsh | Test case archive of optimal power flow (OPF) problems with local optima[END_REF]

Table 1 :

 1 Comparison of bags number I, largest bag size M i and number of variables

	System	I	I by [9] Max. N i	Max. N i by [9] Var. # in (6)	Var. # in (8)	Var. # in (8) by [9]	Found value	Found by [9]
	Polish-2383wp	1242	2383	10	23	2,840,536	23,199	89,893	1.8408 × 10 6	1.8742 × 10 6
	Polish-2736sp	1538	2736	10	23	3,744,216	27,298	104,388	1.3042 × 10 6	1.3082 × 10 6
	Polish-2737sop	1538	2737	10	23	3,746,953	27,034	103,720	7.7572 × 10 5	7.7766 × 10 5
	Polish-2746wop 1546	2746	10	23	3,771,631	29,024	108,950	1.2040 × 10 6	1.2085 × 10 6
	Polish-2746wp	1547	2746	10	24	3,771,631	28,257	107,148	1.6266 × 10 6	1.6324 × 10 6
	Polish-3012wp	1689	3012	10	24	4,537,578	30,996	116,799	2.5727 × 10 6	2.6089 × 10 6
	Polish-3120sp	1757	3120	10	24	4,868,760	32,637	121,869	2.1391 × 10 6	2.1608 × 10 6

Table 2 :

 2 Performance comparison

	System	Found value	Found by [9]	Found by [23]
	Polish-2383wp	1.8408 × 10 6	1.8742 × 10 6	1.8685 × 10 6
	Polish-2736sp	1.3042 × 10 6	1.3082 × 10 6	1.3078 × 10 6
	Polish-2737sop	7.7572 × 10 5	7.7766 × 10 5	7.7763 × 10 5
	Polish-2746wop 1.2040 × 10 6	1.2085 × 10 6	1.2083 × 10 5
	Polish-2746wp	1.6266 × 10 6	1.6324 × 10 6	1.6317 × 10 6
	Polish-3012wp	2.5727 × 10 6	2.6089 × 10 6	2.5917 × 10 6
	Polish-3120sp	2.1391 × 10 6	2.1608 × 10 6	2.1427 × 10 6

Acknowledgements. We thank Cedric Josz, a coauthor of [12] and [18] for sending us these references.


text:
Charles Poussot-Vassal
email: charles.poussot-vassal@onera.fr

David Quero
email: david.queromartin@dlr.de

Pierre Vuillemin
email: pierre.vuillemin@onera.fr

Data-driven approximation of a

Keywords: Model reduction, data-driven, flexible aircraft

published or not. The documents may come

INTRODUCTION

Motivations and contribution 1 : within the aeronautical industry, the design and optimization of an aircraft structure is both a challenging and time consuming task for structural and aeronautical engineers. One of the major specification derived by the authorities before allowing commercial flight is dictated by its response to discrete gust disturbances, modeled by the so-called 1-cosine family gust signals (see e.g. [START_REF] Tang | Experiments and Analysis for a gust generator in a wind tunnel[END_REF]; [START_REF] Haghighat | A model predictive gust load alleviation controller for a highly flexible aircraft[END_REF]). The latter represent a set of specific disturbances exciting the structure around a wide range of frequencies and resulting in an important stress, and varying loads affecting the wingspan. This stress is monitored to control the fatigue and to ensure both passenger safety and aircraft clearance in all flight situations. In addition to the clearance constraint, lightning aircraft structure stands as a major lever to meet aircraft consumption and overall traffic footprint impact reduction objectives. However, the clearance and the lightning objectives are somehow contradictory. Indeed, lightning renders the structure more flexible, leading to different dynamical effects blending2 and might thus render the aircraft more sensitive to gust disturbances (see [START_REF] Poussot-Vassal | Gust load alleviation: identification, control and wind tunnel testing of a 2D aeroelastic airfoil[END_REF]). A challenge for aeronautical, structural and control engineers is then to reduce the overall weight while both minimizing the gust effect, and, in the same time, maintain the aircraft maneuverability (see [START_REF] Meyer | Ground test for vibration control demonstrator[END_REF]). To reach these challenges, a large amount of simulations is usually needed, thus simple but representative models are desirable. Nevertheless, since the involved models are usually highly complex, dedicated numerical software have to be used. Consequently, simulation time and computational burden tend to rise and are becoming a limiting factor. These effects are even more amplified when other engineering fields are involved (e.g. in control, analysis, flight mechanics etc. ) as they may have even higher requirements concerning the complexity of the considered models.

In order to generate a suitable state-space aeroelastic model in the time-domain3 , several approaches have been used. For instance, [START_REF] Smith | Statespace model generation for flexible aircraft[END_REF] provide a summary of the techniques available to generate an aeroelastic model in a state-space form. The most used rational approaches are the ones proposed by [START_REF] Roger | Airplane math modeling methods for active control design[END_REF] and the Minimum-State Method proposed by [START_REF] Karpel | Design for active flutter suppression and gust alleviation using state-space aeroelastic modeling[END_REF]. They both introduce an approximation error due to the basis used to represent the frequency-domain data. The rational approaches have been commonly combined with the Mode Displacement Method (MDM) in order to recover the dynamic loads acting over the airframe [START_REF] Moulin | Gust loads alleviation using special control surfaces[END_REF]. However, the MDM shows a lower convergence rate with an increasing number of structural modes than the Force Summation Method (FSM), see [START_REF] Pototzky | New and existing techniques for dynamic loads analyses of flexible airplanes[END_REF]. In order to avoid the convergence rate problems for the dynamic loads computation, a rational approach can be done on the FSM formulation [START_REF] Pototzky | New and existing techniques for dynamic loads analyses of flexible airplanes[END_REF]). Other authors apply the rational approach over the complete set of panels used for the aerodynamic discretisation (see [START_REF] Castrichini | Gust and ground loads integration for aircraft landing loads prediction[END_REF]), leading to a significant number of additional aerodynamic lag states. Also, additional aspects have to be taken into account when applying the rational approach to a gust input due to exponential complex term in the frequency domain for the associated gust mode. The consideration of several gust modes can avoid such problems at the cost of an increasing number of states, see [START_REF] Wang | Accurate rational function approximation for time-domain gust analysis[END_REF]. Another approach to generate a state-space model of the aeroelastic systems is the Eigensystem Realization Algorithm (ERA), which has been applied to flexible structures (excluding rigid body modes) for the computation of generalized coordinates (see [START_REF] Kim | Aeroelastic model reduction for affordable computational fluid dynamics-based flutter analysis[END_REF]; [START_REF] Silva | Recent enhancements to the development of cfd-based aeroelastic reduced-order models[END_REF]). In this work the ERA has been extended to an aeroelastic system including rigid body modes and dynamic loads as system outputs. However and as stated latter, special care must be taken when applying the algorithm to dynamic loads quantities.

In order to reduce the number of considered states in the timedomain representation of the aeroelastic system, we propose to approximate the HiFi (High Fidelity) model-based simulation results in gust cases (obtained by a complex and accurate simulation setup), by a reduced order LTI (Linear Time Invariant) dynamical model, which is both accurate and simple, and for which extensive simulation can be performed very efficiently. The approach combines techniques from model construction, simulation, and data-driven model approximation with a stable low complexity LTI model. Moreover, the use-case involved in this paper, defined in Section 2, as well as some metrics computation functions will be rendered available as a benchmark for model approximation in a MATLAB format on the MOR Wiki 4 .

Notations: in this paper we denote N, R and C, the natural, real and complex values sets. The H 2 (C + ) (or simply H 2 ) space denotes the set of complex-valued matrix functions H(•) analytic over C + and which inner product integral is bounded along the imaginary axis, the H ∞ (C + ) (or simply H ∞ ) space, denotes the ones which supremum along the imaginary axis is bounded (the H 2 (C -) and H ∞ (C -) spaces are similarly defined). In addition, RH • denotes the spaces of rational functions in H

• . Then L ∞ (ıR) = H ∞ (C + ) ⊕ H ∞ (C -) and L 2 (ıR) = H 2 (C + ) ⊕ H 2 (C -
). H(S)(s) denotes the complexvalued transfer function, where s is the Laplace variable, and S = (E, A, B, C, D) its associated realization of dimension n which eigenvalues λ i , i = 1, . . . , n, are the pencil or generalized eigenvalues of (E, A), and we denote ρ(E, A) = C ∪ {∞} λ i is the resolvant of this matrix pencil. The operator vec stacks the columns of a matrix underneath each other to form a single vector.

Structure of the paper: in Section 2, the aeroelastic structure and aerodynamical dynamic of the considered flexible aircraft benchmark are detailed. The way the data are collected is also described. Then, Section 3 describes the Loewner-based datadriven RH ∞ stable model approximation procedure used for this specific problem, allowing interpolating the data collected with a stable dynamical model. Then, Section 4 illustrates the approach on the considered benchmark and quantify the accuracy control and computational gain. Finally, conclusions and perspectives are discussed in Section 5.

LONGITUDINAL FLEXIBLE AIRCRAFT MODELING

Let us first describe the considered flexible aircraft model and the simulation setup used to obtain the frequency responses, denoted as data later on, on which the model will be constructed. Now an aeroelastic model for perturbations around the nonlinear steady state is obtained. The method is described by [START_REF] Quero | An Aeroelastic Reduced Order Model for Dynamic Response Prediction to Gust Encounters[END_REF] and is based on the unsteady correction of the Aerodynamic Influence Coefficient (AIC) matrices for a set of correction modes. The modified AIC matrices are obtained by solving the linear least squares problem (1),

(W T (k) ⊗ Ã)vec(x) = vec(P Q (M∞, k, αs)) -vec(P Q0 (M∞, k)), (1)
where the symbol ⊗ denotes the Kronecker product and

vec(x) =   Q 11 |Q 11 | • • • Q Nr 1 |Q Nr 1 | • • • Q 1Np Q 1Np • • • Q Nr Np Q Nr Np   T (2)
The number of panels used for the aerodynamic surface discretization is denoted by N p and the reduced frequency k is defined by k = ωL ref /U ∞ , where ω is the circular frequency and L ref the reference length taken as 4.5 (m) for this configuration. The matrix W (k) ∈ C Np×M contains in columns the downwash distribution corresponding to a number M of correction modes, while the matrix P Q (M ∞ , k, α s ) ∈ C Nc×M (where N c represents the number of constrains imposed) contains the set of CFD data for the corresponding correction modes to be matched by the corrected AIC matrices and the matrix A computes the corresponding aerodynamic quantities.

In this application the local lift and local aerodynamic pitch

AIC matrix Q c is Q c (M ∞ , k, α s ) = Q (M ∞ , k) + T Q Q (M ∞ , k, α s ) (3) where matrix T Q properly reshapes the incremental matrix Q from size N r × N p to N p × N p .
Two correction modes (M = 2) corresponding to the gust (assumed to act in the vertical direction and with a constant spanwise intensity) and the first flexible mode are selected for the aeroelastic model generation.The resulting aeroelastic model has been shown to show computational time savings of up to three orders of magnitude when compared to fully coupled CFD-CSM (Computational Fluid Dynamics -Computational Structural Mechanics) simulations. This aeroelastic model produces the HiFi data used in Section 3.

Simulation and data generation: once the aeroelastic model has been applied, a SIMO LTI system with transfer function matrix G can be defined as

G = G T r G T f G T c
T , where

G r = ıωT r H u G f = T f H u G c = ÃT A Q c (D 1jk + ıkD 2jk ) H kg φ gh G u + φ 0 (4) G u = -ω 2 M hh + K hh -q ∞ Q c,hh -1 q ∞ Q c,hj φ 0 φ 0 = n j e z U ∞ e -ıω x j U∞ (5)
where the matrices M hh and K hh denote the generalized mass and stiffness matrices respectively (the structural damping has been neglected), q ∞ the dynamic pressure, D 1jk and D 2jk are substantial derivative matrices, H kg is the spline matrix and φ gh includes the structural modes retained. The matrices Q c , Q c,hh and Q c,hj represent aerodynamic contributions. For a detailed explanation of these matrices see [START_REF] Quero | An Aeroelastic Reduced Order Model for Dynamic Response Prediction to Gust Encounters[END_REF]. The gust mode vector φ 0 ∈ C Np in ( 5) is responsible for the interpolation errors in the classical rational approach due to the presence of the complex exponential term. There, n j is the local normal vector to the panel, e z the (vertical) gust direction and x j the panel x-coordinate defined at the 3/4 local chord point.

The gust input u(ω) is specified by the (nondimensional) gust profile at the aircraft nose and the set of outputs in y

(ω) = [u T h (ω) c T l (ω) c T m (ω)]
T includes the time derivative of the rigid body coordinates, the generalized coordinates corresponding to the flexible modes and the local lift and local pitch aerodynamic coefficients acting over the strips. Note that the time derivative has been considered for the rigid body motion in order to avoid the heave motion to be periodic, which would cause the incremental vertical displacement to return to zero at the end of the simulation. Actually, due to the singularity at zero frequency in the matrix Q c,hh for the column corresponding to the heave motion, the vertical displacement of the aircraft at the end of the simulation is not equal to the initial one for a nonzero gust amplitude. Thus the term ıω in (4) performs the time derivative over the rigid body coordinates. The matrix T r selects the rows corresponding to the rigid body modes and the matrix T f selects the flexible ones. 7) is non unique, then for n, n u , n y ∈ N, by denoting S n,ny,nu (R n×n , R n×n , R n×nu , R ny×n , R ny×nu ), the set of all n y × n u realizations of dimension n, if the matrix pencil (E, A) is regular / non-singular for some finite λ ∈ C,

H(S) : ρ(E, A) → C ny×nu s → C(sE -A) -1 B + D, (8)
is the transfer function associated to S : (E, A, B, C, D) ∈ S n,ny,nu and ρ(E, A) the resolvant of the matrix pencil (E, A).

The set of all realizations of H(S) is denoted S(H(S)). This operator will play an important role in the second part of the paper to ensure model stability.

Preliminaries on data-driven ERA

In the case the E matrix is invertible, the Eigensystem Realization Algorithm (ERA) can also be applied. Thus the first equation in ( 7) can be replaced by ẋ(t) = Āx(t) + Bu(t) after premultiplying by the matrix E -1 . The ERA is then applied over the discrete-time version of the system, x(k + 1) = Ād x(k) + Bd ū(k) and y(k) = Cd x(k) + Dd ū(k), where k = 1, ..., N t with N t the number of time samples. The ERA algorithm is based on the construction of the generalized Hankel matrix containing the discrete-time impulse responses for all input/output combinations and the later singular value decomposition to obtain the state-space representation, see [START_REF] Silva | Recent enhancements to the development of cfd-based aeroelastic reduced-order models[END_REF]. A conversion from the discrete to continuous system allows the final state-space formulation.

In this work the aerodynamic coefficients acting over the component strips are considered in the frequency-domain data, imposing a difficulty in the construction of the impulse responses for the ERA, which is believed to be caused by the addedmass aerodynamic term at increasing frequencies. In order to consider aerodynamic forces in a state-space form where the matrix E is invertible different approaches which explicitly deal with the existence of the added mass term have been presented, see [START_REF] Antoulas | Model reduction and approximation theory and algorithms, chapter A tutorial introduction to the Loewner framework for model reduction[END_REF]. In this work, instead of physically identifying the different contributions, a more general approach is intended by dividing the transfer function data provided in the frequency domain by the term ıω. Consequently, the ū term has to be multiplied by ıω, causing the input to the state-space model obtained by the ERA to be the time derivative of the original input. This introduces a limitation when applying the ERA algorithm to transfer functions including the added-mass term, as the original input variables for controller design and not the time derivatives are sought. Additionally, different subsystem blocks had to be considered in order for the ERA to provide a proper state-space model, leading to an increasing number of states for the free flexible aircraft model. In Section 4, the performance of the (below) proposed method is compared to that obtained with the ERA. Note that the following Loewnerbased approach will allow to get rig of ERA limitations.

Data-driven approximation in the Loewner framework

Let us now assume that data (6) are splitted into two subsets: more specifically, s i into left {µ j } q j=1 and right {λ i } k i=1 and Φ i into left {v j } q j=1 and right {w i } k i=1 . These data are completed with tangential directions {l j } q j=1 and {r i } k i=1 arbitrarily chosen.The Loewner matrices offer a versatile and compliant framework to deal with frequency-domain data, by seeking for a realization interpolating these data, in the barycentric sense. More specifically, let us be given left interpolation driving frequencies {µ j } q j=1 ∈ C with left output or tangential directions {l j } q j=1 ∈ C ny , producing the left responses {v j } q j=1 ∈ C nu and right interpolation driving frequencies {λ i } k i=1 ∈ C with right input or tangential directions {r i } k i=1 ∈ C nu , producing the right responses {w i } k i=1 ∈ C ny , one aims at finding a realization S such that the resulting transfer function H(S) is a tangential interpolant of the data, i.e. satisfies the following left and right interpolation conditions:

l * j H(S)(µ j ) = v * j = l * j Φ j for j = 1, . . . , q , H(S)(λ i )r i = w i = Φ i r i for i = 1, . . . , k . (9)
The main ingredient to achieve ( 9) is the Loewner matrix, which was developed in a series of papers (see e.g. [START_REF] Mayo | A framework for the solution of the generalized realization problem[END_REF]; [START_REF] Antoulas | Model reduction and approximation theory and algorithms, chapter A tutorial introduction to the Loewner framework for model reduction[END_REF]). In the sequel, we recall the main steps. Let be given, the left (or row) data and the right (or column) data:

{µ j , l * j , v * j } q j=1 and {λ i , r i , w i } k i=1 . (10)
Moreover, let us assume that λ i and µ j are distinct, then the associated Loewner L ∈ C q×k and shifted Loewner L σ ∈ C q×k matrices, also referred to as the Loewner pencil, are constructed as follows, for i = 1, . . . , k and j = 1, . . . , q:

[L] j,i = v * j r i -l * j w i µ j -λ i , [L σ ] j,i = µ j v * j r i -λ i l * j w i µ j -λ i . (11)
Then, by organizing the left and right interpolation data as:

M = diag(µ 1 , . . . , µ q ) L * = [l 1 . . . l q ] V * = [v 1 . . . v q ] , Λ = diag(λ 1 , . . . , λ k ) R = [r 1 . . . r k ] W = [w 1 . . . w k ] (12
) the Loewner L and shifted Loewner L σ matrices [START_REF] Silva | Recent enhancements to the development of cfd-based aeroelastic reduced-order models[END_REF] satisfy the following Sylvester equations: LΛ -ML = LW -VR and L σ Λ -ML σ = LWΛ -MVR. Consequently, following the main results of [START_REF] Mayo | A framework for the solution of the generalized realization problem[END_REF], given the right and left interpolation data as in (10), and assuming that k = q and (L, L σ ) is a regular pencil where λ i or µ j are not eigenvalues, the rational transfer function

H(S)(s) = C(sE -A) -1 B, with realization S : (E, A, B, C, 0) constructed as E = -L, A = -L σ , B = V and C = W, (13)
is a minimal descriptor realization which interpolates the left and right constraints, i.e. ensures (9). Now the interpolatory framework, connected to the Loewner one has been reminded, it is obvious that, starting from data (10) either obtained from experiments (see e.g. [START_REF] Poussot-Vassal | Gust load alleviation: identification, control and wind tunnel testing of a 2D aeroelastic airfoil[END_REF]; [START_REF] Meyer | Ground test for vibration control demonstrator[END_REF]) or simulations performed on a complex model, it is simple to construct a rational interpolant S and its associated transfer function H(S) which exactly interpolates the data. One additional possibility embedded in this framework is that it is possible, at reduced numerical cost, to obtain a r-th (r ≤ n) order approximated model. This is obtained through a rank revealing matrix factorisation such as:

L = [ Y 1 Y 2 ] Σ 1 Σ 2 X * 1 X * 2 , (14)
where

Σ 1 ∈ R r×r , Σ 2 ∈ R (n-r)×(n-r) and Y 1 , Y 2 , X 1 , X 2
are of appropriate dimensions. The r-th order reduced order model with realization Ŝ is obtained by Petrov-Galerkin projection:

Ŝ = (-Y * 1 LX 1 , -Y * 1 L σ X 1 , Y * 1 V, WX 1 , 0)
. This later reduce-order model will be exploited in Section 4 using the data collected on the HiFi simulator to derive an approximate model with monitored mismatch error.

Optimal RH ∞ approximation

While the Loewner ensures some interpolatory properties, it does not provide any guarantee concerning the stability of the resulting matrix pencil (E, A). In the perspective of timedomain simulation for aircraft load computations and response to gust disturbance, this is a major issue. Indeed, authorities criterion are based on discrete time-domain gust simulations rather than frequency ones. Therefore, a post treatment should be applied to approximate the (often unstable) Loewner-based interpolant by a stable one. To this aim, the recent work of [START_REF] Kohler | On the closest stable descriptor system in the respective spaces RH 2 and RH ∞[END_REF] is used to project a rational unstable model onto its best stable approximant. Mathematically, let the obtained Loewner-based transfer H(S) ∈ RL ∞ or RL 2 , the objective is to construct H( S) ∈ RH ∞ or RH 2 . In the case where a MIMO LTI continuous-time dynamical system can be represented by a first order descriptor realization S : (E, A, B, C, D) with n u inputs, n y outputs and n internal variables8 , the model is given by a set of DAE as in (7). Similarly to S n,ny,nu , one can define the following sets: (15)

Due to the orthogonality of the L 2 space decomposition from its stable and anti-stable part, the case where p = 2 is straightforwardly solved by simply taking the stable part of H(S). However, in that case, the obtained model is no longer of dimension n, but ñ ≤ n. In [START_REF] Kohler | On the closest stable descriptor system in the respective spaces RH 2 and RH ∞[END_REF], the case where

p = ∞ is solved 9 . By reminding that L ∞ (ıR) = H ∞ (C + ) ⊕ H ∞ (C -),
E = E T -R S-,γ ; A = -A T -R S-,γ -C T -C S-,γ , B = E T -Q -B -; C = C -P -E T -and D = D -, R S-,γ = Q -E -P -E T -γ 2 I (16) evaluated at γ = σ 1 = ||S -|| ∞ ,
and where P -and Q -the solutions of the generalized Lyapunov equations

A -PE T -+ E -PA T -+ B -B T -= 0 A T -QE -+ E T -QA -+ C T -C -= 0. (17)
Moreover, if r = rank(A), there exists the projection matrices P, Q T ∈ R r×n-such that P ÃQ is regular. In this case,

inf G∈ ñ∈N S + ñ,ny ,nu ||H(S) -H(G)|| L∞ = σ 1 . ( 18
)
The proof is detailed in [START_REF] Kohler | On the closest stable descriptor system in the respective spaces RH 2 and RH ∞[END_REF]. The key step consists in performing the stable and unstable part separation and solving two Lyapunov equations. Note that in this case, the resulting stable approximate S and its transfer function H( S) is of dimension ñ = n. Applying this post treatment to the Loewnerbased approximate thus preserves the degree and, hopefully, the accuracy and interpolatory properties.

APPLICATION AND PERFORMANCE EVALUATIONS

Following Section 3, let us now summarize the procedure. It has been clarified that from any frequency-domain measurement set given as (6), one can obtain an exact or approximate interpolant model H(S) with realization S ∈ S 0 n,nu,ny of minimal McMillian degree thanks to the Loewner procedure. Then, in view of time-domain gust simulation, based on the descriptor realization S, the optimal projection onto the RH ∞ space, denoted S ∈ S + n,nu,ny , can be obtained by reasonable algebraic manipulations computations (see Theorem 1), we are now ready to apply it on the considered data generated by the HiFi simulation of the flexible aircraft structure, presented in Section 2. Based on the data collected in Section 2, the Loewner matrices are constructed. Computing the normalized singular values of L, one obtains the singular value decay (not illustrated here space limitations). The rank revealing factorization shows that a perfect match is obtained with n = 321, while more than 420 measurements are used. Constructing H(S) with realization S ∈ S 0 n,nu,ny leads to a perfect matching, but still, with an unusable model (more that 40 unstable eigenvalues). Then, based on the unstable realization, the stability enforcement in the RH ∞ sense, is applied. Then, Figure 3 illustrates the resulting mean and max relative mismatch errors as a function of the approximation order. These errors are the L ∞ and L 2 -norms of the relative mismatch computed between the stable interpolated model H( S)(s i ) for varying order n and the original data Φ i . For n = 100, a mean relative mismatch of 1% is achieved. This error drops around 0.1% for n = 250.

CONCLUSION

In this paper we have presented an aeroelastic benchmark representing a flexible aircraft model which frequency response data have been obtained from a HiFi simulator. Then, as rooted on the Loewner framework, followed with a dedicated postprocessing treatment, ensuring model stability, we obtain a simple but yet accurate model of this complex aircraft usecase, which can be used in place of the original large-scale one for optimisation and intensive simulation. Future works, held in this European project, will include additional inputs of the model (e.g. control surfaces) and a parametrisation of the data as a function of some geometrical data and flight configurations.

4

  http://morwiki.mpi-magdeburg.mpg.de/morwiki/index.php

Fig. 1 .

 1 Fig. 1. M ∞ = 0.84 and α s = 2 (deg). CFD steady pressure distribution (Euler solution). Structural, aerodynamic and aeroelastic models: a configuration based on typical long-range aircraft dimensions is used as reference. The lifting surfaces are defined by NACA0012 profiles. Geometrical and mass properties together with the first natural frequencies of the structure are available in Quero (2017). For the structural model six rigid-body and fifty flexible modes are taken into account. The flexible modes are normalized to yield an identity generalized mass matrix. Note that a very low frequency for the first flexible mode has been considered in order to investigate the applicability of the present approach to flexible configurations. In this work Euler equations discretised over a grid containing approximately 1.5 • 10 5 nodes are considered. The selected transonic flight point is at Mach number 0.84 and a steady angle of attack α s = 2 (deg) corresponding to a trimmed aircraft with total lift coefficient of 0.45 neglecting the effect of the horizontal tail plane on the pitch moment equilibrium. The steady pressure coefficient distribution over the surface as predicted by Computational Fluid Dynamics (CFD) is shown in Figure 1, illustrating a shock wave over the wing component. The corresponding freestream velocity and gust translation speed is U ∞ = 257.93 (m/s).

Fig. 2 .

 2 Fig. 2. Aerodynamic strip over the wing component. moment coefficients acting at the strips of the lifting components have been chosen as shown in Figure 2 and thus the number of constraints N c is twice the number of aerodynamic strips corresponding to the local pitch and local aerodynamic pitch coefficients. Once the unknown vector vec (x) of length N r N p (N r represents the subset of panels which are corrected) is computed, the corrected unsteadyAIC matrix Q c is Q c (M ∞ , k, α s ) = Q (M ∞ , k) + T Q Q (M ∞ , k, α s )(3) where matrix T Q properly reshapes the incremental matrix Q from size N r × N p to N p × N p .

  3. DATA-DRIVEN MODEL APPROXIMATION IN THERH ∞ SPACE3.1 Preliminaries in the data-driven approximationDue to irrational terms in the model description, simulation is performed. Then, based on the frequency-domain collected data, we are now ready to describe the proposed complete twosteps approximation process: the data-driven approximation and projection onto a stable subspace. Within the data-driven approximation framework, we are given{s i , Φ i } N i=1 ,(6)a set of frequency-domain (or complex-domain) n u inputs, n y outputs data Φ i ∈ C ny×nu collected at varying frequencies s i ∈ C either from experimental measurements or any numerical simulation. These data satisfy, for i = 1, . . . , N , y(si ) = G(s i )u(s i ) = Φ i u(s i ),where u(s) ∈ C nu , y(s) ∈ C ny respectively are the Laplace transform values of the inputs u(t) ∈ R nu and outputs y(t) ∈ R ny , evaluated at s i 5 . In this setting, the objective is to find a LTI dynamical model represented by a realization of "complexity" n ∈ N, S : E ẋ(t) = Ax(t) + Bu(t) , y(t) = Cx(t) + Du(t) (7) where x(t) ∈ R n denotes the internal variables 6 and where E, A ∈ R n×n , B ∈ R n×nu , C ∈ R ny×n and D ∈ R ny×nu , are constant matrices, for which the frequency response of the model (7) well reproduces the data set {s i , Φ i } N i=1 7 . Obviously, realization (

  A, B, C, D) ∈ S n,ny,nu | ıR ⊂ ρ(E, A)},S + n,ny,nu {(E, A, B, C, D) ∈ S n,ny,nu | C ≥0 ⊂ ρ(E, A)} and S - n,ny,nu {(E, A, B, C, D) ∈ S n,ny,nu | C ≤0 ⊂ ρ(E, A),E regular}, where S + n,ny,nu and S - n,ny,nu are sets of stable and anti-stable realization-based systems, respectively. Mathematically, given p = {2, ∞} and a realization S ∈ S 0 n,ny,nu i.e. H(S) ∈ RL ∞ or RL 2 , one aims at finding S ∈ S + ñ,ny,nu (ñ ≥ 0) such that, ||H(S) -H( S)|| Hp = inf G∈ ñ∈N S + ñ,ny ,nu ||H(S) -H(G)|| Lp .

  one can similarly write S = S + ⊕ S -, where S ∈ S 0 n,ny,nu ,S + = (E + , A + , B + , C + , D + ) ∈ S + n+,ny,nu and S -= (E -, A -, B -, C -, D -) ∈ S - n-,ny,nu , n = n -+ n + .Then, the following theorem, holds true. Theorem 1. (Optimal RH ∞ approximation). Let be given S ∈ S 0 n,ny,nu (i.e. H(S) ∈ RH ∞ ), S + and S -, then problem (15) for p = ∞ is solved by S = S + ⊕ P S S-,σ 1 Q, where S S-,σ 1 = (E, A, B, C, D) where

Fig. 3 .

 3 Fig. 3. Point-wise relative mismatch error in both the L ∞ (blue rounds) and L 2 (red triangles) sense as a function of the approximation order n.Since the considered data have n u = 1 inputs and n y = 92 outputs, located at different points of the aircraft and wing and representing different measurements (load, moments, etc), it is complicated to illustrate. However, Figure4shows the frequency responses along the right wing span of the HiFi collected data Φ i and the model H( S), for the exact case where n = 321. It clearly shows a very good reproduction / interpolation of the data, by the stable model. Finally, Figure5illustrates the time response of output #1, obtained for varying gust input (left) with the HiFi simulator and the approximated one of dimension n = 100, showing a very good accuracy.

Fig. 4 .Fig. 5 .

 45 Fig. 4. Frequency response along the wing span of the original HiFi model (red dots) and with the exact interpolation model with n = 321 (colored surface).

Usually flight mechanics, loads, vibrations, etc. frequency bands are rather decoupled. When lighter aircraft are designed, these frequency ranges tend to overlap.

Note that time-domain formulation allows for the design of controllers for load alleviation or flutter suppression purposes. However, in this work, the control surfaces are not considered. They can be included in a straightforward way though.

Note that when the data Φ i are collected from experimental or HiFi software simulation, one usually selects s i = ıω i , where ω i ∈ R + is the pulsation of the experiment. The case where s i ∈ C is involved in some cases where the model is analytically known. This is rarely the case when complex simulator enters in the picture.

x(t) ∈ R n are the state variables if E is invertible.

One should note that the above objective is at the boundaries with the traditional model identification problem. However here, an interpolatory framework rather than a norm minimization is sought.

This is the type of model obtained by the Loewner framework.

Note that in both p = 2 and p = ∞ cases one does not consider any eigenvalue on the imaginary axis.

1 This work has been funded within the frame of the Joint Technology Initiative JTI Clean Sky 2, AIRFRAME Integrated Technology Demonstrator platform "AIRFRAME ITD" (contract N CSJU-CS2-GAM-AIR-2014-15-01 Annex 1, Issue B04, October 2nd, 2015) being part of the Horizon 2020 research and Innovation framework programme of the European Commission.