Lucas Richard
email: lucas.richard@g2elab.grenoble-inp.fr

David Frey

Alexis Derbey

Marie-C Écile Alvarez-Herault

Bertrand Raison

Experimental Design of Solar DC Microgrid for the Rural Electrification of Africa

DC microgrids are a promising and cost-effective solution for the rural electrification of Africa or South-East Asia. However, microgrids mainly rely on power electronic converters and those converters must offer, in addition to their normal usage, protection and start-up services as well as low voltage ripple. This must be considered when designing converters for such particular applications. Moreover, field deployment of DC microgrids in rural places is often expensive and logistically difficult. Hence, it is necessary to develop a test bench to test and thoroughly validate the proper operation of a microgrid before its deployment. This paper presents the experimental design needed, from the converter to the test bench, to develop and deploy solar DC microgrids with decentralized production and storage for the rural electrification of Africa.

Introduction

Nowadays, almost one billion people have no or limited access to electricity in the world despite the United Nations (UN) Sustainable Development Goals of ensuring universal access to reliable and modern energy services by 2030 [START_REF]United Nation 17 Goals[END_REF], [2]. The vast majority of those people resides in rural places in Sub-Saharan Africa or South-East Asia, where abundant resources of solar energy are available [START_REF]United Nation 17 Goals[END_REF]- [START_REF] Moner-Girona | Decentralized rural electrification in Kenya: Speeding up universal energy access[END_REF]. Electrifying such rural places through national grid extension is most of the time economically unviable as it requires large upfront costs for cable connections [START_REF] Moner-Girona | Decentralized rural electrification in Kenya: Speeding up universal energy access[END_REF], [START_REF] Nasir | Scalable solar DC microgrids[END_REF]. On the other hand, Solar Home Systems (SHS), i.e. individual solar kits, are undoubtedly gaining momentum in rural places [START_REF] Moner-Girona | Decentralized rural electrification in Kenya: Speeding up universal energy access[END_REF], [START_REF] Groh | Swarm Electrification: Investigating a Paradigm Shift Through the Building of Microgrids Bottom-up[END_REF] but they are only a stopgap measure to quickly give access to basic electricity services to the population. Low-energy access is often combined with high unemployment rate and poor local economic development [START_REF]Nano é presentation website[END_REF], which can not be dealt with SHS. Therefore, to tackle the energy access challenge, it is highly crucial to propose an electrification model which can combine a quick and low-cost access to electricity to inhabitants of rural places and a community uplift through economic development. For instance, Nano é, a French-Malagasy social company created in 2017, proposes an electrification * Institute of Engineering Univ. Grenoble Alpes model based on the construction of electric infrastructures in a bottom-up manner while creating numerous local electric operators [START_REF]Nano é presentation website[END_REF]. Nano é installs and operates solar DC nanogrids (NG), consisting of one solar panel and one lead-acid battery for 4 to 6 houses. Those NGs, a low-cost solution easy to deploy, are only a first step in Nano é's electrification model. In villages where multiple NGs are already installed, it is now highly relevant to interconnect them to form a microgrid (MG) to enhance electricity services, reduce hardware costs and offer productive use of electricity in rural places [START_REF] Nasir | Scalable solar DC microgrids[END_REF], [START_REF] Groh | Swarm Electrification: Investigating a Paradigm Shift Through the Building of Microgrids Bottom-up[END_REF], [START_REF] Nasir | A Decentralized Control Architecture Applied to DC Nanogrid Clusters for Rural Electrification in Developing Regions[END_REF]- [START_REF] Li | A Module-Based Plugn-Play DC Microgrid with Fully Decentralized Control for IEEE Empower a Billion Lives Competition[END_REF]. There is definitely a strong interest for AC and DC solar MGs in the rural electrification sector [START_REF] Nasir | Scalable solar DC microgrids[END_REF], [START_REF] Groh | Swarm Electrification: Investigating a Paradigm Shift Through the Building of Microgrids Bottom-up[END_REF], [START_REF] Li | A Module-Based Plugn-Play DC Microgrid with Fully Decentralized Control for IEEE Empower a Billion Lives Competition[END_REF]- [START_REF] Jhunjhunwala | Solar-DC Microgrid for Indian Homes: A Transforming Power Scenario[END_REF]. However, the topology of such MGs still remains an open question [START_REF] Dragicevic | DC Microgrids -Part II : A Review of Power Architectures[END_REF]- [START_REF] Liptak | Self-organizing nanogrid (SONG)[END_REF]. It is the belief of the authors that a fully decentralized DC topology is preferable, with decentralized production and storage, as advocated by the swarm electrification concept [START_REF] Groh | Swarm Electrification: Investigating a Paradigm Shift Through the Building of Microgrids Bottom-up[END_REF]. Such topology offers a scalable and modular way of building electric infrastructures in a bottom-up manner while being economically viable. In addition, decentralized DC MGs present lower upfront costs and lower losses than centralized topologies [START_REF] Nasir | Scalable solar DC microgrids[END_REF], [START_REF] Nasir | A Decentralized Control Architecture Applied to DC Nanogrid Clusters for Rural Electrification in Developing Regions[END_REF], [START_REF] Nasir | Analysis on central and distributed architectures of solar powered DC microgrids[END_REF]. However, such DC MG entirely relies on power electronic converters and their control. It is therefore needed to develop converters which are suited to those particular applications. There are many types of converters proposed for rural electrification, as the bidirectional buck-boost DC-DC converter [START_REF] Shaw | Modeling and control of a battery connected standalone photovoltaic system[END_REF] or more complex converters (Triple Active Bridge (TAB), flyback) [START_REF] Nasir | A Decentralized Control Architecture Applied to DC Nanogrid Clusters for Rural Electrification in Developing Regions[END_REF], [START_REF] Samende | Power loss analysis of a multiport DC -DC converter for DC grid applications[END_REF]. A converter designed for DC MGs must offer, in addition to its normal usage, protection and start-up services to facilitate the operation of the MGs. Moreover, ripple reduction is particularly important for the stability of the DC bus of the MG as voltage oscillations might interfere with the control algorithms of the MG [START_REF] Dragicevic | DC Microgrids -Part I : A Review of Control Strategies and Stabilization Techniques[END_REF]. An important constraint of rural electrification projects is that field deployment of DC MG pilots in rural areas is often logistically complicated, expensive and time-consuming. Therefore, it is crucial to develop extensive test benches in the lab to thoroughly test the operation of the proposed MGs before any test on the field. This necessitates fast and robust prototyping in the lab to be able to quickly validate a proposed structure without long engineering development time. This paper proposes an experimental design, from the power converter to the test bench, to thoroughly test and validate the proper operation of a DC MG with decentralized storage and production before any field deployment. To the knowledge of the authors, only a few research institutes have designed a test bench for DC MG for rural electrification ( [START_REF] Nasir | A Decentralized Control Architecture Applied to DC Nanogrid Clusters for Rural Electrification in Developing Regions[END_REF], [START_REF] Nasir | Dual-loop control strategy applied to the cluster of multiple nanogrids for rural electrification applications[END_REF], [START_REF] Nasir | Integration and Decentralized Control of Standalone Solar Home Systems for Off-Grid Community Applications[END_REF] at the Center For Research on Microgrids (CROM) facilities and [START_REF] Li | A Module-Based Plugn-Play DC Microgrid with Fully Decentralized Control for IEEE Empower a Billion Lives Competition[END_REF] at the University of Manitoba) and none have had the opportunity to deploy the proposed MG on the field. The rest of the paper is organized as follows. Section 2 presents the design of a converter for a DC MG application. Section 3 details the development of a test bench for the proposed DC MG while section 4 describes the results of the converter and the test bench. Section 5 concludes and quickly introduces the field-test prototype of the converter.

Converter Design for Rural Electrification 2.1 DC Microgrid with Decentralized Production and Storage

The objectives of interconnecting NGs to form a DC MG are manifold. By mutualizing installed production and storage capacities, such a MG would use more efficiently hardware resources, improve the lifetime of the batteries and offer higher electrical services to the community [START_REF] Nasir | Scalable solar DC microgrids[END_REF], [START_REF] Nasir | A Decentralized Control Architecture Applied to DC Nanogrid Clusters for Rural Electrification in Developing Regions[END_REF], [START_REF] Nasir | Dual-loop control strategy applied to the cluster of multiple nanogrids for rural electrification applications[END_REF]. Overall, this would improve the economic viability and sustainability of such projects and help households and communities progressively climb the energy ladder [START_REF] Groh | Swarm Electrification: Investigating a Paradigm Shift Through the Building of Microgrids Bottom-up[END_REF], which would facilitate rapid scaling of those energy-access solutions.

The proposed MG interconnects all NGs in a close area to a 60 V DC bus through bidirectional DC-DC converters, as shown in Fig. 1. Communal loads such as rice huller or water pump can connect to the DC bus and be powered by all NGs. The control algorithm of the converter must be decentralized and communication-free to avoid a single point of failure and be affordably deployable even in areas where telecommunication signals are inexistent or unreliable [START_REF] Nasir | A Decentralized Control Architecture Applied to DC Nanogrid Clusters for Rural Electrification in Developing Regions[END_REF]- [START_REF] Li | A Module-Based Plugn-Play DC Microgrid with Fully Decentralized Control for IEEE Empower a Billion Lives Competition[END_REF]. Moreover, a decentralized and communication-free control algorithm enables plug and play features for the MG [START_REF] Li | A Module-Based Plugn-Play DC Microgrid with Fully Decentralized Control for IEEE Empower a Billion Lives Competition[END_REF], which means that each NG can connect or disconnect to the MG without impacting the operation of the rest of the MG. Therefore, based only on local measurements, the converter controls the power flows between the NG and the MG, i.e. if it must inject or absorb current on the DC bus. A converter assesses the energy state 2 shows the proposed control algorithm, which ensures relevant and consistent power flows on the MG while maintaining the DC bus voltage within a certain range [START_REF] Nasir | A Decentralized Control Architecture Applied to DC Nanogrid Clusters for Rural Electrification in Developing Regions[END_REF], as will be shown in Section 4.2. 3 regions can be defined for the energy state of the NG according to its battery SoC, i.e. low for SoC below 60%, medium for SoC between 60% and 80% and high for SoC above 80%, as defined in Table 1.

Fast Prototyping of a DC-DC Converter

A central element of the DC MG proposed is therefore the bidirectional DC-DC converter interfacing a 12 or 24 V NG to the 60 V DC bus. It has to be highlighted here that rural electrification projects often run on tight budgets and need fast results. Moreover, protection services and start-up procedure are necessary for the proper operation of the MG and it is economically interesting to directly include them in the bidirectional converter. For all those reasons, it has been chosen to design and use a bidirectional buck-boost converter for the proposed MG. Indeed, without galvanic isolation needed, a bidirectional buck-boost converter offers the best compromise between cost, ease and time of design and possible additional services. A bidirectional boost converter would have been sufficient for this application but for start-up procedure and network issues mitigation, a buck-boost topology is preferred, as shown with the schematics of the converter in Fig. 3.

For protection services, mosfets Q1 and Q5 are necessary. Indeed, without mosfets Q1 and Q5, the NG would always be connected to the MG through the diodes of mosfets Q3 and Q7. To properly disconnect the NG to the MG, mosfets Q1 and Q5 must be opened. In addition, Q2 and Q6 are useful to offer a free-wheel path to the current in the inductors L 1 and L 2 in case of over-current and thus limit the constraints on the interface semiconductors (Q1, Q3, Q5, Q7). Moreover, for start-up procedure, the arms Q1/Q2 and Q5/Q6 could be useful. As the DC bus contains high capacitors to stabilize the DC bus voltage and offer some inertia, the initial current at the startup of the first converter can be high to charge all the capacitors to the desired voltage. Therefore, a start-up procedure, where the capacitors on the MG side of the converters are slowly brought to the input voltage thanks to the converter operating in buck mode, might be useful to limit the initial inrush current. Once the output voltage close to the input voltage, the converter can switch to boost operation to bring the DC bus voltage to its normal value.

In addition, under-voltage events on the MG could happen following a fault or a problem with a communal load (high transient current, DC motor suddenly blocked, etc.). If the DC bus voltage suddenly drops to a lower value than the NG voltage, a boost structure would fail due to the presence of the body diode of Q3 and Q7. Such events could therefore not be managed with only a boost converter whereas a buck-boost converter could switch to buck conversion to overcome the under-voltage event or disconnect the NG to the MG by opening Q1 and Q5. Lastly, in anticipation of future usages or features of the MG (e.g. communal loads, higher voltage range of operation), it is of interest to keep the buck operation mode available. By integrating those protection and start-up features, this topology permits to avoid to add costly electromechanical components such as circuit breakers or pre-charge relays to the system which normally perform those functions. Furthermore, as the control algorithm for the converter relies on the DC bus voltage, it is important to obtain low voltage ripple on the DC bus for stability purposes [START_REF] Nasir | A Decentralized Control Architecture Applied to DC Nanogrid Clusters for Rural Electrification in Developing Regions[END_REF], [START_REF] Dragicevic | DC Microgrids -Part I : A Review of Control Strategies and Stabilization Techniques[END_REF]. Therefore, an interleaved converter with 2 arms is chosen to reduce voltage oscillation. Additional CEM filters and a high capacitive bus on the MG side are also included to further reduce voltage oscillations, increase inertia and reinforce stability. For a first prototype, 2 arms were considered to be a good compromise between design simplicity and ripple reduction. Multiple arms also enable modularity and can increase efficiency (only one arm used at low power to decrease switching losses for instance). In addition, it reduces the current per arm for a specified power rating, which eases the design of the power inductors, the rating of the mosfets and the thermal management of the entire converter.

The designed converter can be seen in Fig. 4. As a first step, a prototype mock-up has been conceived in the lab. It is a 18 A converter (9 A per arm), current-controlled with a current sensor in each arm. Its input voltage range is 10 to 29 V (one or two batteries in series), hence it is a 500 W converter. Table 1 gathers information about the proposed converter. The converter is controlled through a myRIO embedded controller from National Instruments programmed in LabVIEW [START_REF]LabVIEW presentation website[END_REF], as it can be seen in

Experimental Test Bench of a DC Microgrid

To validate and further tune the proposed DC MG, it is needed to develop a test bench for thorough testing of the operation of the MG. The objective of the test bench is to create a lab environment as close as possible from the real field conditions, with additional safety, monitoring, ease of installation and testing capacities. The schematic of the test bench is shown in Fig. 5 and the real test bench can be seen in Fig. 6. The test bench represents a MG with 3 NGs connected to each other and consists of 3 bidirectional buck-boost converters connected to a DC bus through RL lines. Each converter is controlled by one myRIO embedded controller. A power supply in parallel to an electronic load can be placed at the input of each converter to emulate a battery with additional features (current and voltage monitoring and current limitation). The RL lines emulate the impedance of the electric cables which would link the NGs, with 4 different available distances (20, 40, 60 and 80 meters) and the cabling set-up can be reconfigured easily as indicated in Fig. 5. A communal load can also be included with a power resistor.

The myRIOs enable to implement many different control features for the DC MG through the LabVIEW software and they are connected to a computer for communication and monitoring purposes. Firstly, the control algorithm of the proposed DC MG (shown in Fig. 2) is implemented on each myRIO. At each converter, the embedded myRIO controller associated measures the current and output voltage of the converter, i.e. the DC bus voltage, then runs the control algorithm which outputs the current reference for the converter. A numeric PI regulator is integrated in the Real-Time module of LabVIEW to control the duty cycle of the mosfets to regulate the current to its reference. Protection features are also included in the myRIOs, which open all high-side mosfets and close all low-side mosfets if any problems occur (such as over-current or overvoltage). In addition, the SoC of the NGs can be emulated on the myRIOs to test different operating points of the MG. The SoC can be automatically changed following a pre-defined scenario or manually changed through a Graphic User Interface (GUI). It is then possible to implement different production/consumption patterns for the NGs, and by using a very small capacity for the emulated batteries, SoC variation due to production/consumption patterns of the NGs can be emulated even on a short timeframe. Lastly, each myRIO records data every 25 ms (voltage, current, current reference, etc.). Therefore, post-treatment of the operation of the test bench is possible, which makes it very insightful. Table 1 summarises the different parameters of the developed test bench. It must be highlighted here that the myRIO embedded controller and labVIEW permit fast prototyping by saving long engineering time associated with the use of micro-controllers and by enabling the implementation of numerous features on the test bench within a relatively small development time. This is believed to be particularly useful for this rural electrification project where the actual goal is to deploy a DC MG on the field and where the test bench is only a intermediary step to validate the proper operation of the proposed MG. Once a proof of concept has been given through a test bench in a lab, it is then relevant to go from a myRIO controller to a micro-controller. 4 Experimental Results

Converter Performances

The converter can operate with one arm or 2 arms interleaved, with the PWM signals of both arms shifted of 180°, as can be seen in Fig. 7. The blue and pink signals show the PWM commands of mosfets Q3 and Q7 whereas the green and yellow signals show the current in power inductors L 1 and L 2 . The currents are shifted by 180°so the output current ripple and then the output voltage ripple are reduced, which is useful for the stability of the control algorithm. It can be noted that the RMS values of both currents are not exactly equal due to small measurement errors which impede the PI regulators of each arm to achieve exact current sharing.

Figure 8 shows the efficiency of the proposed converter for different cases, i.e. one or two arms and 12 or 24 V NG. A higher efficiency is obtained with 24 V NGs than with 12 V NGs, which proves that the main losses are due to conduction losses. For a certain input power, the current will be divided by 2 between a 12 V NG and a 24 V NG which would reduce the conduction losses per a factor around 4 and then increase the efficiency of the converter. As the switching frequency is only 25 kHz, the switching losses are relatively negligible. Also, at a certain power, at a certain voltage (12 or 24 V), there is two times less current in the power inductors with two arms than with one arm. As the conduction losses are quadratic with the current, the efficiency is higher with two arms than with one arm, even at the same voltage, which is confirmed by Fig. 8. With a higher switching frequency and a lower current per arm, this might not be the case and then using only one arm instead of two arms might be relevant to increase efficiency (by turning off one arm, the switching losses of this arm would be 0). Overall, the converter obtains high level of efficiency, especially with two arms with a 24 V NG (up to 99%), which is highly satisfying. If the average operating point of the converter exceeds 150 W, it becomes relevant to change to a 24 V NG to increase the efficiency of the converter. However, it remains an open question whether such a high efficiency traduces or not a global optimum for the design of the converter and the proposed converter does not claim to be perfectly sized. Additional work will be needed in the future to optimise the proposed converter.

On the MG, the first converter to be launched must charge all the DC bus capacitance to the operating voltage. Figure 9 and 10 show the inrush current from the battery and the DC bus voltage when starting a converter respectively in boost mode and in buck then boost mode. This is done with a 90 Ah lead-acid Victron battery charged at 13 V and with only one converter connected to the DC bus. The start-up procedure consists of starting in buck mode to slowly bring the DC bus to a voltage close to the NG voltage, then switch to boost mode. It can be seen that when starting directly in boost mode, the initial inrush current reaches almost 40 A and is quite important for 2 ms whereas when starting in buck mode there is only a 400 mA peak at the beginning in buck mode (as shown by the part zoomed in of Fig. 10) and a 1 A peak at the end of the boost mode, which starts approximately 150 ms after the buck mode. Note that there is a ratio of 40 between the current scales of both figures. The start-up procedure enables to control the current in the buck mode until the DC output voltage reaches close to the NG voltage. Without this start-up procedure, the current can not be controlled until the DC output voltage reaches the input voltage. The initial inrush current without start-up procedure would be even worse with additional capacitance on the DC bus (i.e. with more converters connected to the DC bus) and with a 24 V NG. Furthermore, when a converter starts with an already existing MG (i.e. an other converter is already launched), it is crucial that the converter starts at the proper duty cycle, i.e. V bat /V bus to avoid any high initial current. Different tests have been carried out with a 12 V NG and a 60 V MG (i.e. a proper duty cycle of 20%) and high initial inrush currents have been observed when starting with a duty cycle different from 20%. The converter can even fail and break with an initial duty cycle too far from 20%, which proves that this precaution is essential. Moreover, protection features of the power converter have been thoroughly tested. In case of overcurrent, overvoltage or undervoltage events, the converter reacts by opening Q1 and Q5 and closing Q2 and Q6 (creating a free-wheel path), separating the NG to the MG and discharging the magnetic energy of the inductors within their own resistance. Manual interruption of the converter has also been tested while operating at high power levels, with satisfying results. The protection procedure of opening Q1/Q5 and closing Q2/Q6 enables to limit the constraints on the interface semiconductors, which is believed to preserve the lifetime of those components. The results shown here enable to affirm that the proposed converter presents satisfying performances and that the proposed topology does permit additional protection and start-up features, which is particularly relevant for DC MGs.

Power Flows on the DC Microgrid

The test bench enables to study the operation of the MG. More precisely, the power flows between the different batteries, depending on their SoC level, can be observed as well as the variation of the DC bus voltage. Figure 11 shows the evolution of the DC bus voltage, the emulated SoCs of the 3 batteries and the currents exchanged between the batteries at different operating points. Note that when the current on the NG side is positive, it means that the NG is injecting current on the MG and vice versa. In this scenario, the SoCs are automatically evolving according to different production/consumption patterns implemented in labVIEW. It is important here to remind that a converter determines whether it must inject or absorb current according to the SoC of its NG and the DC bus voltage, as described in Fig. 2. It can be seen that NG 2 starts to inject current to the MG (at t=60 s) as the DC bus voltage is very low and its SoC reaches 60%, and then supports NG 1 and NG 3. As the SoC of the 3 NGs are evolving, the amount of exchanged current reduces to 0 A (t=190 s) when NG 1 and NG 3 reach a SoC of 60%. The DC bus voltage is then stabilized at 60 V which indicates that the global energy level of the MG is medium (i.e., all the batteries have a SoC at a medium level as defined in Section 2.1). Then, due to its production/consumption pattern, NG 2 reaches first a SoC above 80% and starts to inject current to the MG (at t=220 s), supporting NG 1 and 3, until they reach 80% of SoC (at t=470 s). The voltage then stabilizes at 66 V, indicating that the MG is at a high level (i.e., all the batteries have a SoC at a high level). It is interesting to note that when NG 2 starts to inject current at t=220 s, the DC bus voltage shortly stabilizes at 66 V because NG 1 and NG 3 did not respond quickly enough to the increase in voltage. This is due to an inherent delay of the electronic loads used in the test bench when the currents are low. They start to absorb current a few seconds after they should have which makes the DC bus voltage increase momentarily.

Figure 12 shows the operation of the MG with a communal load. The 3 batteries are charged between 60% and 80%, hence the DC bus voltage is at 60 V. At t=10 s, a communal load (i.e. a power resistor in this case) is connected to the DC bus, its power consumption being increased from t=30 s to t=40 s. It is then disconnected at t=200 s and reconnected at t=310 s with a greater power. The 3 batteries inject current in proportion to their SoC level, the greater the SoC, the more current injected. Around t=110 s, the SoCs of the NG 1 and 3 are manually modified to highlight this current sharing feature. Lastly, the DC bus voltage does change with the loading of the communal load, i.e. the greater the power of the communal load, the lower the DC bus voltage. This is logical with the fact that when powering a high power communal load, the MG has in average less energy reserve, hence its DC bus voltage is lower. The operation of the MG with a communal load is satisfying and validated with the results shown in Fig. 12. Many other tests have been successfully carried out (different topologies, meshed, radial, length of lines, etc.) and the results obtained permit to affirm that the proposed DC MG is operating well and that the proposed control algorithm enables consistent power flows while guarantying stability of the MG and maintaining the DC bus voltage within a certain range. It can also be confirmed here that the DC bus voltage does translate in the average energy situation of the MG. This is useful for the use of communal loads or for future interconnections to other MGs or to the national grid [START_REF] Groh | Swarm Electrification: Investigating a Paradigm Shift Through the Building of Microgrids Bottom-up[END_REF], [START_REF] Li | A Module-Based Plugn-Play DC Microgrid with Fully Decentralized Control for IEEE Empower a Billion Lives Competition[END_REF]. For instance, the time at which communal loads are powered could be decided with the DC bus voltage level or the amount of power exchanged between two connected MGs could depend on their respective DC bus voltage level.

Conclusion

This paper presents the development of a DC-DC converter suited for DC MG with decentralized storage and production as well as a test bench of the proposed MG. Performances of the converters and the operation of the test bench are presented in details.

The successful lab prototype of the converter presented here has led to the design of a field-test prototype of an interconnection module, designed to operate autonomously on the field for a long period of time. For practical reasons, additional functionalities have been added to the interconnection module such as LCD screen, SD module, clock, User Interface (UI), as can be seen in Fig. 13 and 14. The interconnection modules have then been installed on the field in Madagascar at the end of 2021 to interconnect 5 NGs and are successfully running since then. This work is only a first step in a design method for power converters suited to their application and future works will focus on thorough study of co-design between the MG and the converter. In particular, as the MG is still under development and open to changes, it is particularly interesting to assess how the MG specifications and the converter design are intertwined to find a technico-economic optimal of the MG/converter system. The second version of the interconnection module will fully take into account those considerations.

Fig. 1 :

 1 Fig. 1: Topology of the proposed microgrid.

Fig. 2 :

 2 Fig. 2: Control algorithm of a converter interfacing a NG to the MG.

Fig. 3 :

 3 Fig. 3: Schematic of the proposed converter.

Fig. 4 :

 4 Fig. 4: Lab prototype of the proposed converter.

Fig. 4 .

 4 The PWM signals driving the mosfets of the converter are generated by the FPGA module of the myRIO, at a frequency of 25 kHz whereas the microprocessor deals with the control algorithm of the converter and the PI regulation.

Fig. 5 :

 5 Fig. 5: Schematic of the proposed test bench.

Fig. 6 :

 6 Fig. 6: Test bench of the proposed microgrid.

Fig. 7 :

 7 Fig. 7: PWM commands and inductor currents of both arms.

Fig. 8 :

 8 Fig. 8: Efficiency of the converter for different cases of usage.

Fig. 9 :

 9 Fig. 9: Initial inrush current without start-up procedure.

Fig. 10 :

 10 Fig. 10: Initial inrush current with start-up procedure.

Fig. 11 :

 11 Fig. 11: Evolution of the voltage, the SoCs and the currents exchanged between the batteries on the test bench.

Fig. 12 :

 12 Fig. 12: Evolution of the voltage, the SoCs and the currents with the introduction of a communal load on the DC bus.

Fig. 13 :

 13 Fig. 13: Command card of the interconnection module.

Fig. 14 :

 14 Fig. 14: Power card of the interconnection module.

Tab. 1:

  Parameters of the converter and the test bench. by measuring the State-of-Charge (SoC) of its battery, and the MG energy state with the DC bus voltage, which represents the global SoC of the MG, i.e the higher the voltage, the more globally charged is the MG and vice versa. High and low limits (±10%) are set on the DC bus voltage. Figure

	Parameters	Symbol Value	Parameters	Symbol	Value
	Number of arms	N	2	Reference voltage of the MG	V ref	60 V
	Power inductor	L 1 , L 2	360 uH Maximum voltage deviation of the MG V min , V max 54 V, 66 V
	DC bus capacitance	C M G	3 mF	High threshold of the battery SoC	SoC max	80%
	Rated current	I rated	18 A	Low threshold of the battery SoC	SoC min	60%
	Input voltage	V bat	10-29 V Lineic resistance of the electric cable R l	1.465 Ω/km
	Power rating	P conv	500 W	Lineic inductance of the electric cable L l	337 µH/km
	Switching frequency	f sw	25 kHz	Proportional and integral parameters	k p , k i	0.005, 0.001
	Number of converters N c	3	PI integration step	δt	200 µs
	of the NG

Acknowledgements

This research, a part of the LEAP-RE project, was funded by the European Union's Horizon 2020 re-search and innovation program under grant agreement N°963530. This paper reflects only the authors' view and the European Climate, Infrastructure and Environment Executive Agency and the European Commission are not responsible for any use that may be made of the information it contains.


domain: spi phys math info

text:
Clément Roos
email: clement.roos@onera.fr

Jean-Marc Biannic
email: jean-marc.biannic@onera.fr

Hélène Evain
email: helene.evain@cnes.fr

A new step towards the integration of probabilistic µ in the aerospace V&V process

Keywords: Probabilistic µ-analysis, Computational tool, Aerospace V&V process

Probabilistic µ-analysis was introduced 20 years ago as a control system validation means able to quantify the probability of rare and potentially critical events. But for a long time, no practical tool offering both good reliability and reasonable computational time was available, making this technique hardly usable in an industrial context. The STOchastic Worst-case Analysis Toolbox (STOWAT) was introduced a few years ago to bridge this gap between theory and practice. It has been significantly improved since then, thanks to the addition of new features, but above all to increasingly efficient implementations, resulting in a dramatic reduction in CPU time. However, until recently, it could only be applied to small-scale models, with up to 4 or 5 uncertainties. In the perspective of analyzing systems with a larger number of uncertain parameters, a time-consuming and tedious process was carried out in the past months. This led to a complete rewrite of the STOWAT, which is now optimized down to the sub-function level, and whose performance is assessed in this paper on a series on benchmarks of increasing complexity with up to about 20 states and 20 uncertainties. This work represents a new step towards the development of a consolidated tool that could reasonably be integrated in the aerospace V&V process in a near future, finding its place between Monte Carlo simulations -useful for quantifying the probability of sufficiently frequent phenomena -and worst-case µ-analysis -relevant for detecting extremely rare events.

Due to their simplicity, MC simulations [START_REF] Helton | Survey of sampling-based methods for uncertainty and sensitivity analysis[END_REF][START_REF] Landau | A guide to Monte Carlo simulations in statistical physics[END_REF] have long been the preferred validation means in the aerospace industry. No analytical representation of the system is needed, and the probability that a whole set of requirements are satisfied is easily computed from time-domain simulations, where the uncertain parameters of the system are sampled based on their probability distributions. But such an approach is generally time-consuming, provide only soft bounds [START_REF] Zhu | Soft vs. hard bounds in probabilistic robustness analysis[END_REF] and may fail in detecting rare events. In contrast with MC simulations, optimization-based techniques intelligently search the parameter space to find an optimum, which makes them suitable for worst-case analysis problems where the aim is to identify rare but possibly critical events. In this context, stochastic optimization methods such as Genetic Algorithms or Cross Entropy based optimization [START_REF] De Boer | A tutorial on the Cross-Entropy method[END_REF] have been shown to perform well in a wide variety of complex aerospace problems [START_REF] Morio | Estimation of rare event probabilities in complex aerospace and other systems: a practical approach[END_REF][START_REF] Wang | Verification and validation framework for autonomous rendezvous systems in terminal phase[END_REF][START_REF] Mujumdar | Cross-entropy based probabilistic analysis of VEGA launcher performance[END_REF], although there are no formal convergence proofs. But their main drawback is the lack of probabilistic bounds, such as the Chernoff bound, to characterize a confidence level.

On the other hand, less expensive deterministic and simulation-free alternatives exist and have reached a good level of maturity, as is the case for µ-analysis [START_REF] Zhou | Robust and optimal control[END_REF][START_REF] Ferreres | A practical approach to robustness analysis with aeronautical applications[END_REF][START_REF] Roos | Systems Modeling, Analysis and Control (SMAC) toolbox: an insight into the robustness analysis library[END_REF]. But unlike MC simulations, if worst-case scenarios are no longer missed, their probability of occurrence is also not measured. Thus, for many problems the worst-case paradigm based on µ can be overly conservative [48] because of the implicit assumption that each of the real uncertain parameters in the model has uniform probability distribution. This can invalidate a control system on the basis of very rare and therefore extremely unlikely events [START_REF] Balas | Analysis of an UAV flight control system using probabilistic µ[END_REF][START_REF] Marcos | Stochastic µ-analysis for launcher thrust vector control systems[END_REF][START_REF] Thai | Probabilistic µ-analysis for stability and H ∞ performance verification[END_REF]. So there is a real need to fill the gap between MC simulations (able to quantify the probability of sufficiently frequent phenomena) and worst-case µ-analysis (relevant to detect extremely rare events), while providing guaranteed bounds on the calculated probabilities which is not possible with Genetic Algorithms or Cross Entropy based optimization. The ambition is to improve the current industrial standard and to fasten the V&V process, which currently accounts for up to 80% of the AOCS total development time in the space industry.

Researchers started in the 1990s to investigate probabilistic µ-analysis [START_REF] Zhu | Soft vs. hard bounds in probabilistic robustness analysis[END_REF][START_REF] Zhu | Improved bounds computation for probabilistic µ[END_REF][START_REF] Khatri | Guaranteed bounds for probabilistic µ[END_REF], which seeks to combine worst-case bounds determined by µ-analysis with probabilistic information. But although attractive from a theoretical perspective, this approach was far from being applicable at that time, as acknowledged by [START_REF] Khatri | Guaranteed bounds for probabilistic µ[END_REF]: It is still not clear how feasible probabilistic µ is. The first dedicated software was developed more than twenty years later [START_REF] Balas | Analysis of an UAV flight control system using probabilistic µ[END_REF]. It was a major improvement, although still not sufficient to address challenging industrial applications, as highlighted in [START_REF] Marcos | Stochastic µ-analysis for launcher thrust vector control systems[END_REF]: In terms of algorithmic implementation, the situation is similar to when first appeared the preliminary implementations of deterministic µ. That is, in some cases it takes prohibitively long and on average takes too long for standard use in the control design cycle -although possibly still acceptable for limited model complexity during analysis in conjunction with Monte Carlo campaigns. The last round began a few years ago, when [START_REF] Falcoz | Probabilistic µ-analysis for system performances assessment[END_REF] developed a new Matlab software based on the SMART Library of the SMAC Toolbox, that implements state-ofthe-art µ-analysis based algorithms to compute robustness margins and performance levels [START_REF] Roos | Systems Modeling, Analysis and Control (SMAC) toolbox: an insight into the robustness analysis library[END_REF]. Shortly afterwards, the first version of the STOchastic Worst-case Analysis Toolbox (STOWAT) was released by [START_REF] Thai | Probabilistic µ-analysis for stability and H ∞ performance verification[END_REF]. Better integrated with the SMART library, it has been significantly improved since then by [START_REF] Biannic | Advanced probabilistic µ-analysis techniques for AOCS validation[END_REF], thanks to the addition of new features, but above all to increasingly efficient implementations, resulting in a dramatic reduction in CPU time. However, until recently, it could only be reasonably applied to small-scale models, with up to 4 or 5 uncertainties. In the perspective of analyzing systems with a larger number of uncertain parameters, a time-consuming and tedious process was carried out in the past months. This led to a complete rewrite of the STOWAT, which is now highly optimized down to the subfunction level, and whose performance is assessed in this paper on a series on benchmarks of increasing complexity. More generally, this work represents a new step towards the development of a consolidated tool that could reasonably be integrated into the aerospace V&V process in a near future. To the best of our knowledge, this is indeed the first time a realistic system with about 20 states and 20 uncertainties -a flexible satellite with sloshing effects in the present case -can be analyzed by a probabilistic µ approach in only a few minutes.

The paper is organized as follows. Section 2 briefly outlines the µ-analysis framework, both from a deterministic and probabilistic point of view. Section 3 presents the latest software tools and shows on a series of benchmarks of increasing complexity how a prototype version limited to academic applications has become in a few years mature enough to integrate an industrial V&V process. These tools are finally applied in Section 4 to the aforementioned satellite benchmark and compared to MC simulations.

µ-analysis framework

Let us consider a continuous-time uncertain LTI system (usually including control laws):

ẋ = A(δ )x + B(δ )u y = C(δ )x + D(δ )u (1)
where the parametric uncertainties δ = (δ 1 , . . . , δ N ) ∈ R N are independent random variables with probability density functions f = ( f 1 , . . . , f N ). It is assumed that A(δ ), B(δ ), C(δ ), D(δ ) are polynomial or rational functions of the δ i . As a result, system (1) can be transformed into a Linear Fractional Representation (LFR) as in Fig. 1 (right): the uncertainties are separated from the nominal (closed-loop) system M(s) and isolated in a block-diagonal operator ∆ = diag(δ 1 I n 1 , . . . , δ N I n N ), where I n i is the n i × n i identity matrix. The set of matrices with the same block-diagonal structure as ∆ is denoted ∆ ∆ ∆. The subset With these notations in mind, two main problems can be solved either in a deterministic framework using classical µ-analysis or in a probabilistic framework using probabilistic µ-analysis. Stability is introduced first. Remark 2 The δ i are often normalized in practice, so that B ∆ ∆ ∆ corresponds to the set of physically meaningful uncertainties. In this case, a natural choice is k = 1 in Problem 1.2 and the associated probability is simply denoted P ∆ ∆ ∆, f (M(s)). This assumption is supposed to be verified in the sequel.

kB ∆ ∆ ∆ ⊂ ∆ ∆ ∆ is then defined as kB ∆ ∆ ∆ = {∆ ∈ ∆ ∆ ∆ : σ (∆) < k} = {∆ ∈ ∆ ∆ ∆ : |δ i | < k, i ∈ [1, N]},
Assuming now that k wc > 1, H ∞ performance is considered next, where T u→y (s, ∆) denotes the transfer from u to y. Once computed, the probabilities P k ∆ ∆ ∆, f (M(s)) and P γ ∆ ∆ ∆, f (M(s)) can be confronted to a given tolerance level ε, so as to validate or reject the considered control system, depending on whether they are lower or higher than ε.

Remark 3

The uncertainties being bounded, their probability distributions must be supported on a bounded interval. Uniform and truncated normal distributions are often used in practice.

The theory behind µ-analysis is not presented in this paper due to space limitations, but the interested reader can for example refer to [START_REF] Zhou | Robust and optimal control[END_REF][START_REF] Ferreres | A practical approach to robustness analysis with aeronautical applications[END_REF][START_REF] Roos | A set of µ-analysis based tools to evaluate the robustness properties of high-dimensional uncertain systems[END_REF] and [START_REF] Thai | Probabilistic µ-analysis for stability and H ∞ performance verification[END_REF][START_REF] Biannic | Advanced probabilistic µ-analysis techniques for AOCS validation[END_REF] for the classical and the probabilistic versions respectively. Only a few facts are briefly recalled below to facilitate the understanding of Sections 3 and 4. First, computing the structured singular value µ is NP-hard in general, so bounds on k wc and γ wc are usually determined in Problems 1.1 and 2.1 instead of the exact values. Much work has been done in the past decades to reduce the gap between these bounds, and (almost) exact values of k wc and γ wc are now obtained in most cases with a reasonable computational time [START_REF] Roos | A detailed comparative analysis of all practical algorithms to compute lower bounds on the structured singular value[END_REF]. Then, probabilistic µ-analysis relies on B&B algorithms to explore the whole uncertainty domain. For stability analysis, this leads to the following partition of the normalized (see Remark 2) uncertainty domain

B δ δ δ = [-1 1] N : B δ δ δ = D s ∪ D s ∪ D s u (2)
where D s , D s and D s u are three sets of disjoint N-cubes corresponding to the domains where stability is guaranteed, instability is guaranteed and stability is undetermined respectively, with probabilities p(D s ), p(D s ) and p(D s u ). The domain D s u stems from the aforementioned NP-hardness issue, but also from the fact that B&B can only approximate D s and D s , and not compute them exactly. The probability p(D s u ) can be reduced by increasing the number of iterations of the algorithm, at the price of an increase in the CPU time. Guaranteed lower and upper bounds on the exact probability P ∆ ∆ ∆, f (M(s)) of instability are finally obtained as follows, thus solving Problem 1.2:

p(D s ) ≤ P ∆ ∆ ∆, f (M(s)) ≤ 1 -p(D s ) = p(D s ) + p(D s u ) (3)
Their accuracy depends on the chosen stopping criterion of the B&B algorithm, which allows to handle the trade-off between accuracy and computational time. Performance analysis can be done in the same way, leading to the following partition:

D s = D γ ∪ D γ ∪ D γ u (4)
where D γ , D γ and D γ u correspond to the domains where performance is guaranteed (

T u→y (s, ∆) ∞ < γ), non-performance is guaranteed ( T u→y (s, ∆) ∞ > γ)
and performance is undetermined respectively, with probabilities p(D γ ), p(D γ ) and p(D γ u ). The following guaranteed bounds on the exact probability

P γ ∆ ∆ ∆, f (M(s)) of non-performance follow: p(D γ ) ≤ P γ ∆ ∆ ∆, f (M(s)) ≤ p(D s ) -p(D γ ) = p(D γ ) + p(D γ u ) (5)
The only difference with stability is that the investigated domain is limited to the domain of guaranteed stability D s , since performance analysis only makes sense for stable systems. The following partition of the normalized uncertainty domain B δ δ δ is finally obtained by combining ( 2) and ( 4):

B δ δ δ = D γ ∪ D γ ∪ D γ u ∪ D s ∪ D s u (6)
3 Computational tools assessment Problems 1.1 and 2.1 can be solved using the SMART Library of the SMAC Toolbox [START_REF] Roos | Systems Modeling, Analysis and Control (SMAC) toolbox: an insight into the robustness analysis library[END_REF][START_REF] Roos | A set of µ-analysis based tools to evaluate the robustness properties of high-dimensional uncertain systems[END_REF], which was introduced in 2013 and implements state-of-the-art µ-analysis based algorithms to compute robustness margins and performance levels. Based on this library, the STOchastic Worst-case Analysis Toolbox (STOWAT) allows to solve Problems 1.2 and 2.2 by computing guaranteed lower and upper bounds on the probabilities P ∆ ∆ ∆, f (M(s)) and P γ ∆ ∆ ∆, f (M(s)) with the desired accuracy. Initially developed under ESA contract RFP/3-16071/19/NL/CRS/hh since 2019, it will be integrated in a forthcoming version of the SMAC Toolbox and made available on the SMAC website w3.onera.fr/smac.

Many improvements have been brought to the STOWAT since the first version presented in [START_REF] Thai | Probabilistic µ-analysis for stability and H ∞ performance verification[END_REF]. But until recently, it could only be reasonably applied to small-scale models. It generally gave good results in the presence of 4 or 5 uncertainties, but beyond that the computation time increased rapidly if a satisfactory accuracy was required. In the perspective of analyzing systems with a larger number of uncertain parameters, a time-consuming and tedious process was carried out in the past months. This led to a complete rewrite of the Matlab code, which is now optimized down to the sub-function level. This tedious work is the price to pay to be competitive with, or even outperform, simulation-based approaches. The result seems to be well worth the effort, as evidenced by the following comparison between the successive versions of the tool.

Let us consider the following simple example extracted from [START_REF] Falcoz | Probabilistic µ-analysis for system performances assessment[END_REF]:

       ẋ = 0 1 -a 1 (δ 1 ) -a 2 (δ 2 ) x + 0 1 u y = 1 0 x (7)
where a 1 (δ 1 ), a 2 (δ 2 ) are two uncertain parameters defined as:

a 1 (δ 1 ) = 1 + 2δ 1 a 2 (δ 2 ) = 0.8 + δ 2 (8)
and δ 1 , δ 2 are two normalized real parametric uncertainties with a uniform distribution on [-1 1]. Four probabilistic µ codes are compared, all implemented with Matlab:

1) the tool described in [START_REF] Falcoz | Probabilistic µ-analysis for system performances assessment[END_REF], which is the first to rely on the SMART Library of the SMAC Toolbox, 2) the prototype version (V0) of the STOWAT introduced in [START_REF] Thai | Probabilistic µ-analysis for stability and H ∞ performance verification[END_REF], which proposes a totally new implementation as well as additional features such as the determination of D s and D γ in addition to D s and D γ , 3) the first consolidated version (V1) of the STOWAT used in [START_REF] Biannic | Advanced probabilistic µ-analysis techniques for AOCS validation[END_REF], directly inspired from the prototype one, but further optimized and robustified, 4) the most recent version (V2) of the STOWAT assessed in the present paper, which consists of a brand new code optimized down to the sub-function level as highlighted above.

Codes #1 and #2 have been kindly provided by D. Alazard and S. Thai respectively, while codes #3 and #4 have been implemented by C. Roos and J-M. Biannic.

The four codes are first applied to solve Problem 1.2 (probabilistic robust stability), and results are presented in Fig. 2 and Table 1, where it can be checked that p(D s ) + p(D s) + p(D s u ) = 100% according to equation [START_REF] Landau | A guide to Monte Carlo simulations in statistical physics[END_REF]. A drastic reduction of CPU time is observed. The curve in Fig. 2 is indeed almost linear on a logarithmic scale and reveals that the computational effort is divided by about 10 from one version to the next. In particular, the CPU time has been divided by:

• 20 between codes #3 and #4, which shows the relevance of the complete rewrite,

• 1200 since the introduction of the STOWAT in 2019, which makes the most recent version now applicable in an industrial context (see also Table 3 and Section 4).

In the meantime, accuracy has also been improved. Code #1 is not able to quantify p(D s) and p(D s u ).

Code #2 can do it, but p(D s u ) remains a bit large for such a simple example. Finally, codes #3 to #4 allow to drastically reduce p(D s u ) from 0.47% to 0.01%, which means that 99.99% of the uncertainty domain can now be categorized in less than 1s. 2, where it can be checked that p(D γ ) + p(D γ ) + p(D γ u ) = p(D s ) according to equation (4). The trend is the same as before, and a drastic reduction of CPU time is observed. In the meantime, accuracy is also improved, with p(D γ u ) decreasing from one version to the next. The decrease in computational time allows to go further in the analysis by choosing a refined stopping condition for the B&B algorithm. This is illustrated in the last column of Table 2, which shows that p(D γ u ) can be divided by almost 6 with a CPU time of only 18s, and in Fig. 3, where the domains of undetermined stability D s u or performance D γ u are hardly visible. is investigated here, but the same trend can be observed for performance analysis. These test cases are characterized by various fields of application, system dimensions and structures of the uncertainties. Some of them contain poorly damped modes, which usually produce extremely sharp peaks on the µ plot, while others have large state vectors as well as numerous and/or highly repeated uncertainties. All results are gathered in Table 3, whose second column deserves some explanation:

• 3/3 means that ∆ is composed of 3 non-repeated uncertain parameters, • 4/5 means that ∆ is composed of 3 non-repeated uncertain parameters and 1 repeated twice, • the drive-by-wire vehicle has 2 non-repeated uncertain parameters and 7 repeated twice,

• the reentry vehicle has 3 non-repeated uncertain parameters and 2 repeated 4 and 6 times, • the hard-disk drive has 19 non-repeated uncertain parameters and 4 repeated twice.

In each case, the stopping criterion is a requested value of p(D s u ) given in column 4, below which the uncertainty domain is no longer divided and the B&B algorithm is interrupted. Note that this threshold is not the smallest that can be reached. It is just selected to have a CPU time of about 10s with code #4. Lower values could probably be reached at the price of an increase of the computational effort, as illustrated for example in Table 2 (see also below).

The results are consistent with those obtained previously on the academic model ( 7)- [START_REF] Zhou | Robust and optimal control[END_REF]. Columns 5 and 6 indeed show that the CPU time is divided by about 10 on average thanks to the complete rewrite. It remains low even for high-order systems, and a good accuracy -measured here by the value of p(D s u ) in column 4 -is obtained in all cases but one. The case of the bank-to-turn missile should be further investigated and reminds us that improvements are probably still possible. Nevertheless, additional tests reveal that values of 6%, 4% and 2% for p(D s u ) can be reached in 11.2s, 24.1s and 112.3s respectively, which remains reasonable. In view of all these results, it appears that the new version of the STOWAT has reached a good level of maturity, which makes it possible to use it on realistic industrial applications, as shown in Section 4. Let us now analyze Table 3 a little further. As expected, good results are obtained even for systems with large state vectors. µ-analysis is indeed a frequency-domain approach, and the size of the frequency response does not depend on the order of the system. It is more surprising that there is no significant increase of the computational time or decrease of the accuracy when the number of uncertainties increases. The most likely explanation is that the B&B algorithm implemented in the STOWAT divides the uncertainty domain along the directions corresponding to the highest µ-sensitivities [START_REF] Braatz | µ-sensitivities as an aid for robust identification[END_REF], i.e. to the uncertainties with the greatest influence on stability. It also shows that the limiting factor of the proposed approach is not the total number of uncertainties as one might think, but the number of uncertainties that most impact stability or performance. It is therefore preferable to have a system with 20 uncertainties, only 4 of which have a high sensitivity, than a system with 8 uncertainties of equivalent sensitivity. Fortunately, it seems that the first case is more frequent than the second one in practice, even if it must be kept in mind that an unfavorable configuration will inevitably occur at some point, which will put the analysis tools in trouble.

Application to a challenging satellite benchmark 4.1 Open-loop satellite modeling

Without a significant loss of generality, since coupling effects are reduced in our context, a singleaxis satellite model is considered next. Let us then denote θ the attitude angle of the satellite about the considered axis and J = 1000 kg.m 2 its main moment of inertia. The rigid dynamic equation simply reduces to:

J θ = T W + T S + T F + T D ( 9
)
where T W is the control input torque (produced by the reaction wheel system), T S is the torque produced by the propellant slosh effects, T F is induced by the flexible modes of the solar arrays attached to the main body and finally, T D captures all remaining disturbance torques (possibly produced by the solar pressure but also by an embedded robotic arm dedicated to on-orbit services). In this work, as is usual in the literature, slosh and flexible effects are represented by poorly damped second-order linear models.

The following generic expression is then obtained:

T S + T F = ∑ i L i s 2 s 2 + 2ξ i ω i s + ω 2 i θ ( 10
)
where the parameters L i , ω i and ξ i respectively denote the magnitude, the frequency and the damping of each mode. The nominal values are presented in Table 4 and the corresponding open-loop Bode plot of the transfer from T W to θ is visualized in Fig. 4.

ξ i 5 × 10 -3 4 × 10 -3 3 × 10 -3 1 × 10 -3 1 × 10 -3
Remark 4 As is further detailed in [START_REF] Biannic | On robust LPV-based observation of fuel slosh dynamics for attitude control design[END_REF], a more realistic control-oriented representation of the slosh effects is obtained with implicitly time-varying frequency and damping characteristics. The latter generally depend indeed on the angular velocity and acceleration of the satellite. Such have been omitted here to generate an uncertain but invariant model as imposed by the µ-analysis framework.

- Finally, as mentioned above, the control input torque T W is generated by a reaction wheel, which can be approximated by a first-order linear model whose uncertain time constant τ ∈ [0.445 , 0.555] will be fixed to its nominal value 0.5 in the design process:

T W = 1 1 + τs T C (11)
Note that this actuator is also rate and magnitude limited. These properties are taken into account in the design phase but not introduced here since the considered maneuvers are assumed to be saturation free.

Structured H ∞ -based robust controller design

A robust attitude tracking controller is designed for the above satellite model. As illustrated in Fig. 5, a second-order reference model R(s) with desired frequency ω r = 0.05 rad/s and damping ξ r = 0.75 is introduced. This model delivers reference signals θ r and θr which are compared to the outputs of the plant. As is clearly visible on the diagram, an extended PID structure is imposed on the controller K(s) which minimizes the H ∞ norm of the weighted transfer from the exogenous inputs w 1 = θ c and w 2 = T D to the exogenous outputs z 1 = ε = θ rθ , its integral z 2 and z 3 = T C . The weighting functions are tuned so that:

• the optimized controller exhibits good robustness properties against actuator delay (delay margin should be above 0.45 s), • the H ∞ norm of the transfer T D (s) = T T D →ε (s) from disturbance torque inputs T D to the attitude error ε should ideally not exceed 0.01.

The last objective which ensures good tracking performance despite perturbations is hard to achieve simultaneously with the delay margin constraint. After a reasonably short tuning process, a fifth-order controller K(s) was yet obtained with the systune routine from the Matlab control toolbox. This controller exhibits a very good delay margin (0.45 s) and acceptably meets the performance requirement with T D (s) ∞ = 1.01 × 10 -2 .

Fig. 5 H ∞ design model

Moreover, parametric robustness properties against possible variations of the inertia, the actuator time constant and the 15 parameters of Table 4 are partly ensured via a multi-model design approach. There is however no guarantee that the whole uncertainty domain is cleared. It is indeed numerically impossible here to include all parametric variations in the design process since there are two many parameters. A validation phase to be performed a posteriori is then required.

Control system validation

Let us now evaluate the pointing performance properties of the closed-loop system by confronting the latest version of the STOWAT (code #4 in Section 3) to a classical MC-based approach. Real parametric uncertainties are introduced on the inertia J of the satellite and the time constant τ of the reaction wheel, as well as the magnitudes L i , damping factors ξ i and frequencies ω i of the 5 slosh and flexible modes. An uncertain system with 21 states and 17 uncertainties is therefore obtained. Problem 1.1 is solved first using the SMART Library, which shows that stability is guaranteed when all uncertainties are allowed to vary independently of each other by ±11.34% around their nominal values. The main objective is now to determine if the pointing performance degradation with respect to the nominal case is acceptable when the uncertain parameters vary by ±11%. For this, the transfer T D (s, ∆) between the external disturbance torque T D and the attitude error ε is considered. The nominal H ∞ norm without uncertainties is 1.01 × 10 -2 , which indicates a good rejection of the perturbation. Robust performance is first analyzed using random sampling. 100000 samples are generated for each uncertain parameter, following a uniform distribution on an interval of ±11% around their nominal value. The highest H ∞ norm among these 100000 samples is 1.11 × 10 -2 and is achieved near the frequency ω s = 0.5 rad/s. This suggests that performance is only minimally affected by the uncertainties. The Bode plots of the nominal system and the 100 samples with the highest H ∞ norm are shown in red and blue respectively in Fig. 6. There is a fairly large dispersion on the frequency, but much less on the magnitude, which is consistent with the previous results. All computations are performed in 1300 s. The SMART Library is now used to solve Problem 2.1, i.e. to compute a lower bound on the worst-case H ∞ norm of the transfer T D (s, ∆). 100 initial frequencies are considered to cover the whole frequency range, and are then optimized during the analysis process to find the highest possible lower bounds corresponding to worst-case uncertainty combinations. The results are obtained after only 9 s and are partially reproduced in Table 5. It can be seen that in most cases, a lower bound of about 1.1 × 10 -2 is obtained close to ω s , which is in line with the results obtained previously by sampling. But iteration 93 reveals a much higher bound of 2.89 × 10 -2 at the frequency ω c = 2.346 rad/s, which is missed by the sampling approach despite the fact that 100000 cases are considered. The Bode plot of this worstcase configuration is plotted in yellow in Fig. 6. It is worth noting that ω c is quite different from ω s . Moreover, it does not coincide with any frequency of the slosh and flexible modes (0.1, 0.2, 0.3, 0.6 and 1 rad/s). It was therefore not possible before applying µ-analysis to anticipate this phenomenon, which can have a significant impact on the pointing performance. Let us indeed consider a sinusoidal disturbance T D (t) = sin (w c t). The resulting pointing error is plotted in Fig. 7 with the same color code as above. For the worst-case configuration identified by the SMART Library, it is about 15 times and 30 times larger than for the worst random sample and the nominal configuration respectively. In this context, a critical issue is to compute the probability of observing such a performance degradation, so as to decide whether the proposed control system is reliable enough or not. The latest version of the STOWAT introduced in Section 3 is therefore applied. As with the previous sampling approach, all uncertainties follow a uniform distribution. But the frequency range on which performance is investigated is now limited to the interval [2 3] rad/s (represented by vertical dashed lines in Fig. 6) to focus on the worst-case peak identified above. Problem 2.2 is solved with a desired performance level γ = 1.01 × 10 -2 , which corresponds to the H ∞ norm of the nominal system. In other words, the objective is to compute guaranteed probabilities that the H ∞ norm is lower or higher between 2 and 3 rad/s than its nominal value on the whole frequency range (represented by a horizontal green line in Fig. 6). The probability of undetermined performance below which the B&B algorithm should be interrupted is fixed to p(D γ u ) = 0.001%. After about 100 s, the probabilities p(D γ ) = 99.999%, p(D γ ) = 0% and p(D γ u ) = 0.001% of guaranteed performance, guaranteed non-performance and undetermined performance are obtained, which means that the probability that the H ∞ norm is higher than 1.01 × 10 -2 between 2 and 3 rad/s is upper bounded by 0.001%. A very large part of the uncertainty domain has therefore been cleared in a very reasonable computational time. The results are even better if we consider a more realistic case where the uniform distribution is replaced with a normal distribution with zero mean and variance σ 2 . To illustrate this, all probabilities are quickly recomputed using the STOWAT for different values of σ 2 between 0.1 and 5, without performing the whole analysis again. The results are plotted in Fig. 8 and show that p(D γ u ) decreases when σ 2 decreases, the uniform case being recovered when σ → ∞. For example, if σ 2 = 0.2, the probability of non-performance is no larger than 1.8 × 10 -6 %. And results are even better if σ 2 = 0.1, since the probability of non-performance is bounded by 1.3 × 10 -9 %, which is negligible. The performance degradation shown in Fig. 6 and7 is therefore very unlikely in these cases. To sum up, if the system is well enough identified to ensure that σ remains quite small, there is certainly no need to invalidate the control system on the basis of this extremely rare worst-case event. But if σ is larger, performance degradation cannot be ignored, and either a better controller should be designed or the system should be better identified. This realistic example clearly shows the added value of the probabilistic µ approach implemented in the STOWAT, which is here faster and more reliable than MC simulations, but also less conservative than deterministic µ-analysis. Moreover, the computational time remains particularly low -a few minutes to build Fig. 8 -considering the complexity of the model, which contains almost 20 uncertainties.

Conclusion

The combination of state-of-the-art algorithms and a highly optimized implementation reveals the interesting potential of probabilistic µ to evaluate the robustness properties of complex systems. The latest version of the Matlab STOchastic Worst-case Analysis Toolbox indeed allowed to assess the pointing performance of a flexible satellite subject to slosh effects in the presence of about twenty parametric uncertainties, whereas previously existing tools could hardly handle more than 4 or 5. The implementation effort reported in this paper paves the way to the integration of probabilistic µ in the aerospace V&V process in complement to traditional MC-based techniques.

While it is always possible to further optimize the code, the effort tends to become inversely proportional to the benefit, and it is questionable whether it is worth it. More promising directions to make the analysis faster would be to use parallel computing, and to improve the B&B strategies to explore the uncertainty domain more efficiently. Research is currently underway and suggests that by combining the two, it should be possible to further reduce the computational time by a significant amount, leading to a third generation of the STOWAT capable of handling even more complex systems.

  SMAC = Systems Modeling Analysis and Control SMART = Skew Mu Analysis Robustness Tools STOWAT = STOchastic Worst-case Analysis Toolbox V&V = Verification and Validation 1 Introduction

Remark 1 Fig. 1

 11 Fig. 1 Standard interconnections for robust stability (left) and performance (right) analysis

Problem 1 . 1 (

 11 Deterministic worst-case stability) Compute the largest value k wc such that the interconnection of Fig. 1 (left) is stable for all ∆ ∈ k wc B ∆ ∆ ∆ . Problem 1.2 (Probabilistic robust stability) Given a desired stability margin k ≥ k wc , compute the probability P k ∆ ∆ ∆, f (M(s)) that the interconnection of Fig. 1 (left) is unstable when ∆ ∈ kB ∆ ∆ ∆ .

Problem 2 . 1 (

 21 Deterministic worst-case H ∞ performance) Compute the smallest value γ wc > 0 such that T u→y (s, ∆) ∞ < γ wc on Fig. 1 (right) for all ∆ ∈ B ∆ ∆ ∆ . Problem 2.2 (Probabilistic robust H ∞ performance) Given a desired performance level γ ∈ [0 γ wc ], compute the probability P γ ∆ ∆ ∆, f (M(s)) that T u→y (s, ∆) ∞ > γ on Fig. 1 (right) when ∆ ∈ B ∆ ∆ ∆ .

Fig. 2

 2 Fig. 2 CPU time for the successive versions of the tool (stability)

Fig. 3

 3 Fig. 3 Partition of the uncertainty domain with a refined stopping criterion

Fig. 4

 4 Fig. 4 Bode plot of the nominal open-loop transfer function of the satellite model from the control input torque T W to the attitude angle θ (with 3 slosh modes and 2 flexible modes)

Fig. 6

 6 Fig.6Bode plot of the transfer T D (s, ∆) between the external disturbance torque and the attitude error (red = nominal, blue = 100 worst samples, yellow = worst-case identified by µ-analysis)

Fig. 7

 7 Fig. 7 Pointing error in the presence of a sinusoidal disturbance T D (t) = sin (w c t) (red = nominal, blue = 100 worst samples, yellow = worst-case identified by µ-analysis)

Fig. 8

 8 Fig. 8 Probability of undetermined performance as a function of the variance σ 2

Table 1

 1 Numerical results for the successive versions of the tool (stability)

	Code #	1	2	3	4
	CPU time [s] 1208	144	20	<1
	p(D s ) [%]	92.52 92.59 92.79 92.79
	p(D s) [%]	-	6.94	7.20	7.20
	p(D s u ) [%]	-	0.47	0.01	0.01
	The four codes are applied next to solve Problem 2.2 (probabilistic robust performance), and results
	are summarized in Table

Table 2

 2 Numerical results for the successive versions of the tool (performance)

	Code #	1	2	3	4
	CPU time [s] 3370	132	27	<2	18
	p(D γ ) [%]	41.30 41.30 41.33 41.57 42.30
	p(D γ ) [%]	-	48.68 49.40 49.39 50.18
	p(D γ u ) [%]	-	2.61	2.06	1.82	0.32

Codes #3 (V1) and #4 (V2) of the STOWAT are now applied to a series of study cases of increasing complexity extracted from the SMAC Toolbox

[START_REF] Roos | A detailed comparative analysis of all practical algorithms to compute lower bounds on the structured singular value[END_REF]

. Only probabilistic stability analysis (Problem 1.2)

Table 3

 3 Application of the STOWAT to study cases of increasing complexity

	Benchmark	Nb of uncertainties Nb of Requested CPU time [s] CPU time [s] and size of ∆ states p(D s u )[%] Code #3 Code #4
	Inverted pendulum	3/3	4	0.8	32.3	1.9
	Anti-aliasing filter	4/5	2	0.05	116.7	9.9
	Bank-to-turn missile	4/4	6	8	102.8	7.4
	Cassini spacecraft	4/4	17	0.06	35.1	3.2
	Drive-by-wire vehicle	9/16	4	0	17.5	3.9
	Reentry vehicle	5/13	7	0	19.7	2.5
	Space shuttle	9/9	34	0.5	50.7	5.3
	Hard-disk drive	23/27	29	1	102.3	19.1

Table 4

 4 Slosh & flexible modes coefficients (nominal values)

	Mode #i	1	2	3	4	5
	L i [kg.m 2 ]	30	40	50	300	100
	ω i [rad/s]	0.1	0.2	0.3	0.6	1

Table 5

 5 Worst-case performance analysis results Grid point Initial frequency [rad/s] Final frequency [rad/s] Lower bound on T D (s, ∆) ∞

	1	0.000	0.475	1.10 × 10 -2
	2	0.004	0.483	1.10 × 10 -2
	91	1.772	0.489	1.10 × 10 -2
	92	2.070	0.484	1.10 × 10 -2
	93	2.418	2.346	2.89 × 10 -2
	94	2.825	0.476	1.10 × 10 -2
	99	6.141	0.491	1.11 × 10 -2
	100	7.173	0.499	1.11 × 10 -2