Sergey Klyapovskiy

Shi You

Andrea Michiorri

George Kariniotakis

Henrik W Bindner

Incorporating flexibility options into distribution grid reinforcement planning: A techno-economic framework approach

Keywords: Active distribution network, flexibility characterization framework, active elements, flexibility services. Nomenclature

Distributed energy resources (DER) and new types of consumer equipment create many challenges for the distribution system operators (DSOs). Power congestions that can potentially be created during normal or contingency situations will lead to increased investments into grid reinforcement. An alternative solution is to use the flexibility provided by the local resources in the grid. In this paper value of flexibility (VoF) is used as an indicator that can be utilized by the DSO to compare it against costs of the different active elements (AEs) providing flexibility services (FSs). The paper proposes flexibility characterization framework that allows to generalize the process of the cost estimations of any AE by using combinations of cost functions. A case study based on an actual distribution grid is provided to demonstrate

q d dielectric loss, [W/m] r resistance of the conductor, [Ω/m] R T,1 , R T,2 , R T,3 thermal resistance of cable's insulating layers, [m • C/W ] R T soil thermal resistance, [m • C/W ] T temperature ratio, [1/ • C]
T cable , T cable,rated , T cable,i,mean cable's current, rated and mean temperature,

[ • C]
T internal , T external temperatures of internal and external parts of the cable,

[ • C] T S soil temperature, [ • C] z soil layer, [cm]

Introduction

The intensifying efforts to reduce CO 2 emissions and increase energy efficiency are bringing large number of distributed energy resources (DER), such as distributed generation (DG) [START_REF] Jenkins | Distributed Generation[END_REF][START_REF] Strachan | Distributed generation and distribution utilities[END_REF], energy storage systems (ESS) [START_REF] Saboori | Energy storage planning in electric power distribution networks-a state-of-the-art review[END_REF][START_REF] Chouhan | Review of energy storage systems[END_REF][START_REF] Geth | Integration of energy storage in distribution grids[END_REF] and new types of consumer devices like heat pumps (HPs) [START_REF] Akmal | Impact of heat pump load on distribution networks[END_REF] and electric vehicles (EVs) [START_REF] Haidar | Technical challenges for electric power industries due to grid-integrated electric vehicles in low volt-age distributions: A review[END_REF][START_REF] Zhao | Maximizing the penetration of plug-in electric vehicles in distribution network[END_REF] to the distribution grids. The presence of these new components is changing the way distribution networks are operated. The electric power congestions [START_REF] Yusoff | Congestion management in power system: A review[END_REF][START_REF] Haque | Congestion management in smart distribution network[END_REF] or out-of-range voltages [START_REF] Zhao | Understanding lv network voltage distribution-uk smart grid demonstration experience[END_REF] that could be predicted in the previously passive networks now appear with a very short warning, requiring immediate actions.

In order to make the distribution networks more flexible and resilient towards demand profiles with high uncertainties and rapid variations, one of the most prominent solution for the distribution system operator (DSO) is to use flexibility [START_REF] Holttinen | The flexibility workout: managing variable resources and assessing the need for power system modification[END_REF] provided by the consumers or grid's components, known as active elements (AEs) [START_REF] Spiliotis | Demand flexibility versus physical network expansions in distribution grids[END_REF][START_REF] Klyapovskiy | Utilizing flexibility services from a large heat pump to postpone grid reinforcement[END_REF]. Devices like HPs, EVs, ESS or even circuit-breakers (CBs) and transformers with on-load tap changer (OLTC) can be treated as AEs. The services they offer to the system are flexibility services (FSs) [START_REF] Mathieu | Coordination of flexibility services within a distribution network[END_REF], since they allow the distribution system to better adapt to the current power demand situation and potentially evolve into an active distribution network (ADN) [START_REF] Adamo | Global survey on planning and operationof active distribution networks update of cigre c6.11 working groupactivities[END_REF]. Long-term planning of the distribution networks, which before was considered almost independent of the system's daily operation, will become deeply intertwined with the operational planning, if the DSO chooses to utilize FSs [START_REF] Liu | Review and prospect of active distribution system planning[END_REF][START_REF] Pilo | New electricity distribution network planning approaches for integrating renewable[END_REF][START_REF] Kazmi | Smart distribution networks: A review of modern distribution concepts from a planning perspective[END_REF].

Different AEs possess different characteristics and can be operated in various ways, while still providing the same FSs, like voltage regulation or peak reduction. Therefore distribution grid planners could potentially face the situations, where there is a need to choose between different AEs to achieve the same planning objective. In order to do that, AEs should be compared on the same basis. From the DSO perspective, such common ground in most of the cases is the cost of implementing the solution.

In most of the cases, FSs from AEs, that are used to solve distribution network problems are alternatives to conventional solutions, such as grid reinforcement, load curtailment, installation of voltage correction equipment, etc. The value of flexibility (VoF) is the maximum cost that the DSO is ready to pay for the FSs and is determined for each case individually. If the cost of using AEs exceeds the VoF, contracting the corresponding FSs will not be economically viable for the DSO.

The paper proposes a flexibility characterization framework that provides a generic approach to the cost estimation of typical AEs, while taking into account AEs' distinctive features. Quantification of the benefits and expenses from each AE through the framework will enable the comparison between different AEs and conventional solution and help the distribution planner in the related decision-making process.

The paper is structured as follows: Section II contains the state-of-theart concerning typical grid issues, quantification of the FSs, description of AEs used in the paper and their cost estimation. Four types of AEs/FSs are considered: demand response (DR), battery energy storage system (BESS), reconfiguration (RE) and dynamic line rating (DLR). Section III describes the methodology for determining the VoF. The general description of the flexibility characterization framework is given in Section IV. Section V presents the case study highlighting the potential application of the proposed framework, the results of which are shown in Section VI. Finally, conclusions are drawn in Section VII.

State-of-the-art

Distribution grid issues

AEs could provide flexibility, that can be used by the DSO to solve issues in the distribution grid.

The most common issues at the distribution level are power congestions and out-of-range voltages [START_REF] Karagiannopoulos | Hybrid approach for planning and operating active distribution grids[END_REF][START_REF] Huang | Review of congestion management methods for distribution networks with high penetration of distributed energy resources[END_REF]. The out-of-range voltages could be either under-or overvoltages caused by various reasons, such as voltage drop, DGs (e.g. PV panels), motor start-ups, etc. Power congestion occur when the power flowing through a grid component exceeds its maximum carrying capacity (MCC), thus blocking a part of the network. Power congestion issues will be used as an example throughout this paper.

Typically, distribution grid components, such as transformers and cables allow short-term overloading [START_REF] Nawzad | Short-time overloading of power transformers[END_REF][START_REF] Swingler | Statistics of ac underground cables in power networks[END_REF] without jeopardizing the system performance, which is used by the DSOs. That is why from the DSO's experience congestions are more likely to occur during contingency situations, such as N-1 [START_REF] Xiao | N-1 loadability for distribution systems[END_REF]. N-1 contingency is a system state, where one of the grid's components is out of operation due to a fault or a failure. To minimize the loss of load in N-1, some of the remaining grid components will have to accommodate the increased power flows. With the increasing number of consumer devices with the cyclic operation modes (HPs, electric heaters, EVs, etc.), it could be expected to have power congestions even in the normal operation, if a large number of such devices are to turn on simultaneously.

Planning horizons used by most of the DSOs lie within 10-20 yearslong-term planning. However, estimation of the power consumption and parameters of the potential congestion events (CEs) requires high accuracy forecasting, which is problematic to achieve within a long time range. To be able to utilize the benefits provided by the FSs of AEs more efficiently, the planning horizons could potentially be reduced to a medium-term range (ca.

5 years).

Quantification of FSs

Before a DSO can use FSs from AEs in planning, the amounts of available FSs should be properly quantified. In [START_REF] Sun | Benefits of smart control of hybrid heat pumps: An analysis of field trial data[END_REF] the data-driven approach is suggested to estimate the amount of FSs obtained from a HP. By analyzing the data from the two test trials, the amounts of FSs provided in response to different tariff schemes are mapped and generalized with results applicable to the other HP systems. Similar approaches based on the data from smart grid projects are presented in [START_REF] Bernards | Metaanalysis of the results of european smart grid projects to quantify residential flexibility[END_REF][START_REF] Ginsberg | Grid flexibility and the cost of integrating variable renewable energy: Toward a renewable energy integration adder for san diego gas and electric service territory and the california electric grid[END_REF]. The dynamic flexibility function is discussed in [START_REF] Junker | Characterizing the energy flexibility of buildings and districts[END_REF]. The flexibility function shows the response of a smart building to a penalty signal like price or CO 2 level. The characteristics of such a function could be used for the flexibility quantification.

The probabilistic methodology involving clustering techniques is described in [START_REF] Kouzelis | Probabilistic quantification of potentially flexible residential demand[END_REF]. First, all consumers are divided into flexible and non-flexible clusters. Then the probability distribution of non-flexible consumers is built. By comparing it to the power consumption patterns of flexible consumers, the amounts of potential flexibility could be quantified. Optimization algorithms for node flexibility estimation are shown in [START_REF] Ji | Quantified flexibility evaluation of soft open points to improve distributed generator penetration in active distribution networks based on difference-of-convex programming[END_REF][START_REF] Ji | Quantified analysis method for operational flexibility of active distribution networks with high penetration of distributed generators[END_REF]. The node flexibility is defined as a relation between the generation and consumption at each node.

The objective function attempts to balance each node by providing flexibility from different devices without violating their operational constraints. The methodology could be applied to increase the penetration of DER. Finally, [START_REF] Oluleye | A two-step optimization model for quantifying the flexibility potential of power-toheat systems in dwellings[END_REF] proposes the two-step optimization model to quantify the flexibility of power-to-heat systems. The flexibility is quantified based on how much extra power production could be absorbed by the system.

Cost estimation of AEs

In order to obtain the most cost-efficient planning solution, DSO has to compare different planning alternatives such as traditional reinforcement and using FSs from AEs. A comprehensive review of planning algorithms involving AEs could be found in [START_REF] Klyapovskiy | Incorporate flexibility in distribution grid planning through a framework solution[END_REF].

Cost estimation of AEs is based on an understanding of what factors form

the final cost and allow the DSO to decide whether it is economically viable to use AEs in planning. The literature review of the cost models for DR, BESS, RE and DLR is given below.

DR

Consumer devices are typically turned on based on the user's necessity or convenience. In the case of the latter, it is possible to postpone the moment when the device starts to operate. The change in the device's operational pattern caused by an external signal from the DSO or an aggregator is called DR. DR is an FS that can be provided by a wide range of electric devices (EVs, HPs, electric heaters, electric boilers, lighting, etc.). [START_REF] Zvirgzdina | Aggregator as cost optimization tool for energy demand[END_REF] shows the DR cost model based on the supply-demand curves made by an aggregator. [START_REF] Gutiérrez-Alcaraz | Effects of demand response programs on distribution system operation[END_REF] describes how the capacity of the DR FSs is changing following the corresponding change in price signal. The change depends on the self-and cross-elasticity coefficients of each type of DR. Calculating the bid price of the DR by taking users' comfort level into account is described in [START_REF] Chassin | A new thermostat for real-time price demand response: Cost, comfort and energy impacts of discrete-time control without deadband[END_REF]. In [START_REF] Abbasi | Risk-constrained offering strategies for a price-maker demand response aggregator[END_REF] stochastic mix-integer optimization algorithm is proposed to optimize the cost of acquiring DR FSs for an aggregator. [START_REF] Domínguez | The economic impact of demand response on distribution network planning[END_REF][START_REF] Li | Dispatch model of active distribution network based on demand response[END_REF][START_REF] Gu | Dynamic pricing for responsive demand to increase distribution network efficiency[END_REF] show the cost models based on the fixed incentives, average or dynamic pricing. More examples of the methodologies for DR modelling are given in [START_REF] Klaassen | A methodology to assess demand response benefits from a system perspective: A dutch case study[END_REF][START_REF] Boudec | Demand response using service curves[END_REF][START_REF] Reihani | A novel approach using flexible scheduling and aggregation to optimize demand response in the developing interactive grid market architecture[END_REF].

BESS

BESS is a device for storing electrical energy, which can be used at a later point in time. From the DSO's perspective, BESS could be considered as a local generator, able to provide power to the nearby electrical loads during the CEs in the distribution networks, thus potentially alleviating components' overloading. BESS is an AE, typically owned by an independent party (e.g.

aggregator) due to the unbundling rules in the power systems [START_REF]Status Review on the Implementation of Distribution System Operators Unbundling Provisions ofthe 3rdEnergy Package[END_REF].

[45] describes the BESS cost model based on the levelized cost of storage.

The model estimates both capital (CAPEX) and operational (OPEX) expenditures of using BESS. The cost of recharging cycles that are diminishing the BESS resource is included in the OPEX estimation. In [START_REF]Battery energy storage system (bess): A cost/benefit analysis for a pv power station[END_REF] the detailed analysis of the BESS CAPEX that consists of storage unit cost, power conversion system and inefficiency factor is given. Inefficiency factor represents a relation between the rated and actual energy that could be extracted from the BESS. OPEX model of the BESS is shown in [START_REF] Ahsan | Dynamic battery operational cost modeling for energy dispatch[END_REF]. The model is split into an electrical part, where the state of charge is determined and a degradation part which estimates the reduction of a unit's lifetime (LT). Other cost models for the BESS are described in [START_REF] Nguyen | Stochastic optimization of renewable-based microgrid operation incorporating battery operating cost[END_REF][START_REF] Zhou | Modeling of the cost of ev battery wear due to v2g application in power systems[END_REF][START_REF] Yang | Cost-benefit study of dispersed battery storage to increase penetration of photovoltaic systems on distribution feeders[END_REF][START_REF] Correa | Optimal scheduling of storage devices in smart buildings including battery cycling[END_REF][START_REF] Marchi | Life cycle cost analysis for bess optimal sizing[END_REF][START_REF] Choi | Bess life span evaluation in terms of battery wear through operation examples of bess for frequency regulation[END_REF].

RE

Process of changing system topology via altering its power supply routes is called RE. By using RE in the distribution grid, it is possible to shift the electrical load from the congested part of the network to the part with the spare capacity, thus ensuring that all electrical customers are provided with power, and system itself is not in jeopardy. RE is a FS provided by the CBs. Since CBs are an integral part of the distribution networks, CBs are the example of utility-owned AEs.

[54] determines the cost of RE via the cost of switchings (SWs). The model is further expanded in [START_REF] Liu | A network reconfiguration method considering data uncertainties in smart distribution networks[END_REF] by adding the CB's operation time to the number of SWs. [START_REF] Milani | A new probabilistic approach for distribution network reconfiguration: Applicability to real networks[END_REF][START_REF] Tzeng | Generic switching actions of distribution system operation using dynamic programming method[END_REF][START_REF] Ghasemi | Radial distribution systems reconfiguration considering power losses cost and damage cost due to power supply interruption of consumers[END_REF][START_REF] Vaskantiras | Value assessment of distribution network reconfiguration: A danish case study[END_REF] present similar cost models for utilizing RE in planning.

DLR

In a traditional operation of distribution networks, the power ratings of overhead lines and cables are considered as fixed values independent of the external conditions. However, the exact MCCs of the lines are constantly changing and determined by the ability to dissipate the heat created by an electric current. Dynamically changing MCC based on the external conditions is referred to as DLR. By applying DLR in the ADN, a larger degree of overloading could be allowed, thus potentially eliminating CEs. DLR is a FS provided by the lines (overhead or cables), which are utility-owned AEs.

[60] describes the methodology for identifying the benefits and costs of using DLR. The cost is the cost of sensors (e.g. for soil temperature) needed to enable DLR, while the benefits are the amounts of extra energy that could flow through the cable. The DLR cost model presented in [START_REF] Teng | Understanding the benefits of dynamic line rating under multiple sources of uncertainty[END_REF] shows how adding the forecast of DLR capabilities to the real-time measurements might increase the value of the DLR FSs. The degradation model of the cables under different loading conditions is given in [START_REF] El Aziz | Estimation of the lifetime of electrical components in distribution networks[END_REF]. More DLR cost models could be found in [START_REF] Michiorri | Investigation into the influence of environmental conditions on power system ratings[END_REF][START_REF] Bracale | A new procedure to forecast the dynamic rating of buried cables[END_REF][START_REF] Van Genuchten | A closed-form equation for predicting the hydraulic conductivity of unsaturated soils 1[END_REF][START_REF] Buhari | Modelling of ageing distribution cable for replacement planning[END_REF].

The majority of the existing literature regarding the cost estimation and quantification of AEs is focusing on a specific technology under the specific conditions. To facilitate the integration of AEs in planning, DSO needs a generic framework, that could be used to characterize the cost of any AE using the set of the same criteria.

Value of flexibility

VoF represents the highest price DSO is ready to pay to the FS providers for their aid in solving a specific issue in a distribution grid. Since VoF depends on the forecasting of the power demand, in the current paper a probabilistic approach is applied to evaluate its value under different demand forecasting scenarios in the long-term planning. The percentile P % will produce the value, which is higher than or equal to the values in the sample with P % probability [START_REF]Percentiles, percentile rank & percentile range: Definition & examples[END_REF]. By choosing the 95th percentile, the values for the worst case could be evaluated (highest cost, largest power demand, etc.).

VoF can be determined as the minimum cost of the DSO's conventional planning alternatives. In case of potential congestion in the distribution network, DSO can either reinforce the grid or curtail the loads. In the latter case, DSO will have to pay the penalty for the energy not supplied (ENS), but it will also get the benefits from the saved investments in the reinforcement (e.g. interest in the bank). VoF can be calculated as follows:

V oF y,P % = min

     C conv,T OT EX,y,P % , C EN S * E EN S,y,P % -BF y,P % , (1)
Since the reinforcement of components will allow the planners to solve the problem in the long-term perspective (assuming the forecast and planning solution are adequate), to correctly compare flexibility and conventional planning, the annuity is used to calculate the cost of the conventional solution as will be shown in Eq. 8. The annual cost of the conventional solution depends on the LT of the components that have to be installed. Depending on their expected LT (e.g. LT of some cables could be 40 years), the cost of reinforced components per year in most of the situations would be lower than the cost of ENS. VoF for the next year will depend on the decisions made in the previous years. By reinforcing the grid, both the cost of ENS and the cost of conventional solution in the next year will be lower, thus reducing the VoF for that year.

The benefits from postponing the reinforcements by doing load curtailment BF y,P % in Eq. 1 are estimated using the cost of conventional solution without annuities (i.e. total cost). This is done in order to evaluate the amounts of income that could be generated, if the whole amount of C conv,T OT EX,y,P % would be placed in the bank with certain interest rate:

BF y,P % = C conv,T OT EX,y,P % * dr (2)
As could be seen from Eq. 1, TOTEX C conv,T OT EX of the conventional planning solution and the amounts of E EN S are required to calculate VoF.

By comparing VoF with the cost of flexibility planning C f lex,T OT EX , the most cost-efficient solution is the one providing the highest benefits according to Eq. 3:

BF f lex,y,P % = V oF y,P % -C f lex,T OT EX,y,P %

Benefits in Eq. 3 could take positive or negative values. Positive ben- Both conventional and flexibility planning types are described below.

Conventional planning

To illustrate how the cost of a conventional planning solution is obtained,

Fig. 1 showing a potential congestion situation during N-1 is considered.

Based on the forecasts made for each year until the end of the planning horizon pl, in some of the scenarios SC power flow on one of the feeders in the distribution system is expected to exceed cable's MCC S cable,M CC during N-1 contingency with estimated peak demand:

S est,max,peak,N 1,SC = max 1≤y≤pl S est,max,N 1,y,SC , (4)
The conventional solution, in this case, will be to reinforce that feeder with a cable with higher MCC according to the dimensioning criteria S dc,conv , which is determined using the following expression:

S dc,conv ≥ 1 L * S est,max,peak,N 1,P % , (5)
Parameter L links component's rated power and its MCC. Short-term overloading coefficient is either obtained from the component's specification or determined by the DSO based on its own experience. Different values of L could be used for different durations of the overload, the shorter it is, the higher the value of L. For the long-term planning purposes, L corresponds to the duration of the longest possible overloading event (i.e. in the range of hours).

TOTEX of any planning solution for a given year consists of two parts -CAPEX and OPEX according to the following equation:

C T OT EX,y,P % = C CAP EX,y,P % + C OP EX,y,P % , (6)
C CAP EX,y,P % and C OP EX,y,P % are calculated using C CAP EX,y,SC and C OP EX,y,SC for all planning scenarios and value of percentile P %.

If CAPEX and OPEX occurring at different years have to be summed together, the net present value (NPV) formula should be used. NPV allows to refer all the expenses to a current year:

C SC = pl y=0 C y,SC (1 + dr) y , (7)
The annuity payments from the fixed costs of the equipment, that are dependent on the LT of components, are used to calculate C conv,CAP EX,y,SC in Eq. 8: To get the cost of conventional solution without annuities C conv,T OT EX,y,SC for Eq. 2, the term LT comp,n in Eq. 8 should be omitted.

C conv,CAP EX,y,SC = f (S dc,conv ) = (8) = N C n=1 C comp,
OPEX are expenses that depend on the number of hours components are in operation and are distributed throughout the components' LT. OPEX will vary for each scenario. Future OPEX for the year y for scenario SC:

C conv,OP EX,y,SC = f (S est,y,SC , h comp,n,y,SC ) = (9) = N C n=1 C comp,M &L,n,y,SC ,
As was shown in Eq. 1, DSO could choose to curtail the loads and accept the penalty for E EN S,y,P % . In this paper E EN S,y,P % is defined as the area between the S cable,M CC and S est,max,peak,N 1,P % , i.e. between the current network capacity and forecasted maximum power demand (Fig. 1). E EN S,y,P % is equal to the E f lex,y,P % , that is described in the next section.

Flexibility planning

An alternative to the conventional approach is flexibility planning. It utilizes the potential of FSs from AEs to achieve peak reduction during the congestion hours as shown in Fig. 2. The dimensioning criteria for the flexibility planning will be S dc,f lex :

S dc,f lex ≥ 1 L * S target , (10)
The power demand exceeding S target should be covered by the FSs. Eq. 10 is given for the case, when only one AE is used to provide all the flexibility needed. In the situations when multiple AEs are used, additional subindex k denoting the kth AE should be used together with updated power demand curves S est,k,y,SC that takes into account the effect of each AE and different values of S target for each k, so that:

If S target is chosen to be equal to S cable,M CC , S dc,f lex will be equal to S cable,rated and cable reinforcement could be postponed.

In the flexibility planning, the periods when CEs are supposed to occur (shaded areas on Fig. 2) at each year y of the planning horizon should be analysed and the following parameters determined for each scenario SC: duration of using flexibility during each CE D f lex,i,y,SC , capacity S f lex,i,y,SC

and energy E f lex,i,y,SC requested from the flexibility for each CE i and total number of times flexibility should be used N O f lex,y,SC during year y. These parameters are used for determining the VoF and the costs of using AEs.

The calculation of C f lex,CAP EX,y,SC and C f lex,OP EX,y,SC parts of the flexibility solution is described in the next section.

Flexibility characterization framework

Proposed flexibility characterization framework is shown in Fig. 3 C1 identifies whether an AE has to be built prior to the provision of FSs, parameter C1 = 1 if construction is required and 0 otherwise. Based on the answer to the question G2 about AE's ownership, this cost may be included in the CAPEX;

Additional equipment required? (C2)

Additional equipment such as ICT infrastructure for DR or transformer substation expansion to fit the transformer with an OLTC could be required to enable the provision of the FSs from a specific AE. Parameter C2 = 1 if there is a need for additional equipment and 0 otherwise.

OPEX parameters

OPEX part of the total cost of the flexibility planning C f lex,T OT EX,y,P % is determined through "OPEX parameters". It is proposed to estimate OPEX via cost functions corresponding to the questions below. The cost from each cost function is determined for each planning scenario, the final value is then chosen using the P % percentile.

How much and how long? (O1)

The cost of the requested energy from FS for ith CE can be described by O4 determines how much in advance warning has to be given to the AE for providing FS: F O4 (N P f lex,k,i,y,SC ).

a

Evaluation criteria

Total cost of a flexibility planning solution with one or several AEs C f lex,T OT EX,y,P % for a year is used as an evaluation criteria (TOTEX of AE (E)) in the proposed framework. It consists of the CAPEX and OPEX parts, that are described below.

Based on the proposed flexibility characterization framework, CAPEX of the flexibility solution C f lex,CAP EX,y,SC can be calculated as follows:

C f lex,CAP EX,y,SC = K SC k=1, G1 k ,G3 k =0 C f lex,CAP EX,k,y,SC = (12) K SC k=1, G1 k ,G3 k =0 (G2 k * C1 k * C AE,k,y,SC + C2 k * A AE,k,y,SC ),
C f lex,CAP EX,y,P % is calculated using C f lex,CAP EX,y,SC for all planning scenarios and the P % percentile (Eq. 6-7). The cost C AE,k,y,SC is a function of the chosen S target .

OPEX from all used AEs C f lex,OP EX,y,SC at the year y in scenario SC depends on the following cost functions and parameters:

C f lex,OP EX,y,SC = K SC k=1 C f lex,OP EX,k,y,SC = (13) = K SC k=1 f (G1 k , G2 k , G3 k , C(F O1 k,y ) SC , C(F O2 k,y ) SC , C(F O3 k,y ) SC , C(F O4 k,y ) SC ),
In the current paper, the following expression is used to represent the relationship between all the parameters and cost functions for OPEX calculation:

C f lex,OP EX,y,SC = K SC k=1, G1 k ,G3 k =0 C f lex,OP EX,k,y,SC = (14)
=

K SC k=1, G1 k ,G3 k =0 [C(F O1 k,y ) SC + C(F O2 k,y ) SC + + C(F O3 k,y ) SC + C(F O4 k,y ) SC + C AE,f eat,k,y,SC + + G2 k * C AE,M &L,k,y,SC ],
Similar to the CAPEX part, C f lex,OP EX,y,P % is calculated using C f lex,OP EX,y,SC and Eq. 6-7.

Table 1 shows what costs are included in CAPEX and OPEX for each of the four considered AEs/FSs. CAPEX costs are related to the installations of the main and additional components, where the main component is an AE itself. OPEX would typically include the cost of the supplied energy or the cost of the extra degradation caused by using a particular AE/FS. The detailed explanation of how the cost of one of the FS -DLR is modelled is given in the Appendix as an example.

For BESS, RE and DLR the costs could be to a certain degree linked to the physical processes occurring in the devices, while the cost of the DR is more subjective. An extra cost that should be added to the DR is the difference in the energy cost at original time (when the customer intends to use equipment) and later time (when equipment is actually used due to postponement). This would ensure, that the DR owner will not have to pay higher energy costs for postponing its power consumption. An extra cost for the BESS includes a cost of a number of recharging cycles used to provide requested energy (since each cycle will cause BESS degradation) and the partial cost recovery of a BESS unit cost.

Case study

System topology

The system considered in a case study is shown in Fig. 4. It is a part of the real 10 kV distribution network of Nordhavn area in Copenhagen, Denmark. The system is supplied from the 30/10 kV main substation (MS) through the four main cables MS-1, MS-10, MS-2, and MS-20. The topology is organized in two loops with the possibility to shift the electrical load from one feeder to the other in case of a failure or a fault (internal RE via TL T19-T110 or T27-T28).

Each bus represents a 10 kV side of the secondary substation with one or several distribution transformers. The loads are connected to the 0.4 kV side of the transformers (Fig. 5). The power demand at each substation 4 is formed by a combination of various residential, commercial and light industrial consumers. The time-series consumption data with 1-h resolution is synthesized using the annual energy measurements from actual Nordhavn consumers for the year 2016 and corresponding demand curves for each load category. Maximum power of the substations is given in Table AI in the Appendix.

Flexibility sources

The system has a number of FSs providers. DR and BESS utilize the consumers' equipment present in the grid to provide FSs. RE and DLR use the distribution network components as the source of FSs.

Loads at the 0.4 kV nodes 111L1, 114L1, 24L1, 26L1, 29L1, 211L1 and 214L1 can provide DR. In the current paper the capacity of DR is considered as a percentage of the total substation load and change from year to year with the upper limit of 10% of the total load. BESS units of various capacity and energy are connected through their own designated transformers to the 10 kV nodes of 15, 113, 26, 212 and 213.

Parameters of BESS are given in Table AII in Appendix.

Flexibility could be obtained by shifting part of the electrical loads to the feeders in another loop (external RE). To enable external RE three TLs and six CBs have to be constructed (shown in green in Fig. 4).

By installing sensors to measure soil temperature and making a thermal model of the cable it is possible to apply DLR at any cable. The extra capacity given by the DLR is changing from year to year depending on the temperature and the number of rainfall with the upper limit of 15% for the cables in clay.

In order to facilitate the integration of the FSs providers in the longterm planning, without compromising DSO's performance in a benchmarking [START_REF]Incentives Schemes for Regulating Distribution System Operators[END_REF][START_REF]The challenges of grid loss in the regulation of danish distribution system operators. danish energy regulatory authority[END_REF] and its ability to supply customers with power, it is assumed that the FSs from AEs are provided based on the long-term contracts. Those contracts specify the maximum amount of provided capacity, the maximum duration of FSs and limit the number of times, services are provided during a specified time interval (e.g. a year). It is assumed that there is no market, where the DSO can contract the FSs from.

Results

Time-series load profiles for the year 2016 are used as a base for creating forecasts for the period of 4 years with an yearly interval. Each load is decomposed into a base trend, seasonal variation and stochastic component.

At each year ten different scenarios of how the base trend and the seasonal variations will change are considered, resulting in a total of 100 forecasts.

Example of the different scenarios of power demand evolution is given in Fig.

A1 in Appendix.

The analysis of the forecasted scenarios shows that cables C2-21, C214- 213 and C20-214 could experience overload in year 3. Cable C213-212 could also be overloaded in year 4 in addition to the previous three. The overload could potentially occur during the worst case of N-1 contingency (fault at one of the main cables in a loop, when the remaining main cable takes the whole load of the loop). The information about detected CEs is summarized in Table 2. 34, respectively.

The entire process of changing the cables is assumed to take one year, therefore the investment decision should be taken one year before the CEs could occur. The cost of the cables in Table 3 include the cost of components and their installation.

The costs C conv,CAP EX,y,95 , C conv,OP EX,y,95 and C conv,T OT EX,y,95 shown in Table 4 are the reinforcement costs for the cables. Three cables: C2-21, C214-213 and C20-214 should be reinforced in year 3 and one cable C213-212 -in year 4. VoF for the year 3 is calculated using Eq. 1 assuming all three cables would be reinforced in the conventional planning solution. As was mentioned earlier, VoF in year 4 will depend on the decisions made in year 3. Therefore, two values for the VoF are given for year 4. The first value assumes that the DSO decides to reinforce the cables in year 3 (cables C2-21, C214-213 and C20-214). In that case, the cost of the conventional solution The second value of VoF is calculated for the case, when the DSO decides to postpone the reinforcement of cables in year 3. The cost of conventional planning solution now is the sum of reinforcement costs of cables in year 3 and 4. As could be seen from Table 4, the longer reinforcement is deferred the higher is the VoF.

Flexibility planning

Since the probability that the worst case of N-1 will coincide with the high power demand is not high, it may be more economically beneficial to rely on the flexibility options that require much lower CAPEX investments.

The flexibility planning is used to identify which AEs could potentially help in handling detected CEs. The calculation of the total cost of AE is done as described in the flexibility characterization framework using cost functions as the one for RE in Fig. 6. Other cost functions for DR, BESS and DLR could be found in Appendix (Fig. A2 -A8) with parameters used to define them given in Table 5.

The maximum capacity and energy that could be provided by each AE and the costs of utilizing that AE are summarized in invested in the reinforcements in year 3, the VoF in year 4 is only 0.8 ke and therefore using AEs will bring the benefits of 0.5 ke (saving 62.5% of the expenses in conventional planning); if the DSO did not upgrade its network in year 3, the VoF will become higher and therefore the benefits will reach 15.2 ke (saving 98.0% of the expenses in case of reinforcement). As could be seen, involving the FSs from AEs in the distribution network planning can make the network more cost-efficient, while keeping its performance on an acceptable level.

Conclusions and future work

The paper proposes the flexibility characterization framework. The aim of the framework is to generalize the process of cost estimation of various AEs/FSs that could be used as an alternative to the traditional grid reinforcement. Comparing TOTEX from using AEs with the VoF, the DSO planners will have a useful tool to aid them in the decision making process.

Provided case study demonstrates the application of the framework in the part of an actual 10 kV distribution system of Nordhavn area in Copenhagen.

Proposed framework could be used by the DSOs to decide which AEs to choose for providing FSs in their networks. The methodology for the cost estimation of the AEs and VoF could be applied by the DSOs to design the price signals for their DR programs. Potential AEs such as electric heat boosters tested in EnergyLab Nordhavn project [START_REF] Cai | Technical assessment of electric heat boosters in low-temperature district heating based on combined heat and power analysis[END_REF] could be integrated into distribution network planning with the use of the presented framework. By streamlining the process of integrating new flexibility that can provide load reduction as well as load increase, a larger number of renewables could be connected to the electrical network.

In regards to the future work, the identification of the parameters used in the cost functions could be improved by analyzing data from the DSOs that are already using different types of AEs in practice. In this paper, the methodology for VoF estimation is given from the perspective of the DSO.

However, the VoF is different for different stakeholders involved in the network operation (e.g. aggregator). This should be further investigated. The DSO is not the only entity that might require FSs from AEs in its network.

In the future, a competition for the FSs between a DSO and transmission system operator could be expected, which will influence the prices the DSO should offer to the owners of AEs.

  f lex benefits of postponing the reinforcement and benefits of O1 ), C(F O2 ), C(F O3 ), C(F O4 ) OPEX from using AE, determined by the cost functions, [e] C1, C2 CAPEX parameters of the framework, [-] C AE,f eat , C AE,M &L additional cost due to the unique features of an AE and cost of maintenance and losses caused by an AE, [e] C AE , A AE self-cost and additional expenses caused by installation of an AE, [e] C comp,M &L cost of the maintenance and losses for the component, [e] C comp , IC comp , A comp self-cost, installation cost and additional expenses caused by the installation of the component, [e] C conv,CAP EX , C conv,OP EX , C conv,T OT EX CAPEX, OPEX, TOTEX of the conventional planning solution, [e] C conv,T OT EX TOTEX of the conventional planning solution without annuity, [e] C EN S cost per unit of ENS, [e/kVAh] C f lex,CAP EX , C f lex,OP EX , C f lex,T OT EX CAPEX, OPEX, TOTEX of the flexibility planning solution, [e]

Figure 1 :

 1 Figure 1: Estimated worst case power demand on the feeder at year y for scenario SCconventional planning to solve power congestion on a feeder

Figure 2 :

 2 Figure 2: Estimated worst case power demand on the feeder at year y for scenario SCflexibility planning to solve power congestion on a feeder. CE -congestion event

Figure 4 :Figure 5 :

 45 Figure 4: Part of the 10 kV distribution system of Nordhavn area used in the case study, nodes with loads are 0.4 kV. Figure explanations: MS -main substation; L -load node; T -node with TL; CB -circuit-breaker; nodes with the red font have BESS; nodes with the orange font have DR; new supply paths to enable RE with an external loop are highlighted in green

Figure 6 :

 6 Figure 6: Cost function FO2 for RE in years 3 and 4

  S est,max,N 1,y , S est,max,peak,N 1 estimated peak power demand during N-1 at the year y and during the whole planning horizon pl,[kVA]S est estimated worst case power demand for specific operation mode (e.g.

	D f lex duration of using an AE, ∆T difference of temperatures, nc number of conductors in the cable,	[h] [ • C] [-]
	dr δ θ , δ T soil hydraulic and soil thermal diffusivity, discount or interest rate,	[-] [m 2 /s]
	E EN S total amount of ENS, normal operation, N-1, etc.), λ 1 , λ 2 ratios between losses in the metal sheaths of the cable and its total [kVAh] [kVAh] [kVA] losses, [-] E f lex total requested energy from the flexibility, G1, G2, G3 general parameters of the framework, [-] S f lex , S f lex,mean capacity and mean capacity requested from an AE during CE, [kVA] ρ T soil thermal resistivity, [m • C/W ]
	h comp number of hours the component is in operation, S target target peak power limit maintained throughout the planning horizon, [h] σ S,dry dry-soil density, [kg/m 3 ]
	K θ	total number of AEs used for providing flexibility, moisture content,	[-] [kVA] [-]
	L Subscripts short-term overloading coefficient, B degradation process coefficient,	[-] [ • C]
	LT comp component's LT, i ith CE C cable,total total cost of the cable,	[e] [e]
	N C total number of components, k kth AE C LT cable cost of an hour of cable's LT,	[-] [e/h]
	N O f lex total number of times flexibility from an AE is used, etc. n nth individual distribution network component: cable, transformer, [-] C ST soil specific heat, [W s/kg • C]
	N P f lex notification period for requesting FS from an AE, P % percentile cosφ power factor,	[h] [-]
	pl SC SCth scenario planning horizon, I cable DLR of the cable,	[year] [A]
	ROC f lex rate-of-change of requested capacity from an AE, y yth year j number of conductors in the cable,	[kVA/h] [-]
	S cable,M CC cable's MCC, Appendix parameters k θ soil hydraulic conductivity,	[kVA] [m/s]
	S cable,rated cable's rated apparent power, α attainment factor, LT cable,D f lex,i,start cable's LT at the start of ith CE with T cable,i,mean ,	[-] [kVA] [h]
	C T OT EX TOTEX of a planning solution, ∆LT cable , ∆LT cable reduction in cable's LT due to a CE with T cable,i,mean and [e] S dc,conv , S dc,f lex dimensioning criteria for the component in the conventional LT cable,rated cable's LT at rated conditions, [h]
	D f lex,end ending point of AE's usage during CE, and flexibility planning, reduction referred to T cable,rated , N soil composition,	[h] [kVA] [h] [-]

  n,y,SC +IC comp,n,y,SC +A comp,n,y,SC LTcomp,n

, Additional expenses A comp,n,y,SC from the installation of the nth component could include expanding the substation in case of a new transformer, new ICT network, etc.

  cost function F O1 (S f lex,k,i,y,mean,SC , D f lex,k,i,y,SC );

	4.3.2. How often? (O2)
	F

O2 (N O f lex,k,y,SC ) is a cost function showing the cost of total number of times FS is used; 4.3.3. How fast? (O3) Cost function F O3 (ROC f lex,k,i,y,SC ) determines the cost of AE's ramping up/down capabilities in providing FSs; 4.3.4. How long in advance? (O4)

Table 1 :

 1 Costs included in the CAPEX and OPEX of different AE/FSs according to the

	flexibility characterization framework
			Costs included in CAPEX	Costs included in OPEX
	AE/FS	AE					Cost functions
			Main	Additional	Main	Additional
	DR	Consumer's	-	ICT and control	Supplied energy Energy difference	O1-O4
		equipment		infrastructure
						Recharge cycles,
	BESS	BESS	-	ICT and control	Supplied energy	recovery cost	O1-O2
				infrastructure
						of BESS unit
	RE	CB	CBs	Tie-lines (TLs)	Number of SWs	-	O2
	DLR	Cable	-	Soil temperature	Reduction in LT	-	O1
				sensors

Table 2 :

 2 Detected CEs in cables. Values for power, energy, number of CEs and duration are obtained by using 95th percentile

			Year 3					Year 4
		SC						SC
	Cable	with	Sflex,3,95,	Eflex,3,95,	N Oflex,3,95,	Dflex,3,95,	Cable	with	Sflex,4,95,	Eflex,4,95,	N Oflex,4,95,	Dflex,4,95,
			[kVA]	[kVAh]	[-]	[h]			[kVA]	[kVAh]	[-]	[h]
		CE, [-]						CE, [-]
	C2-21	99	1334.4	3578.5	8	2	C2-21	99	1538.3	4867.9	11	2
	C214-213	93	975.8	1543.9	4	2	C213-212	16	166.7	166.7	1	1
	C20-214	99	1334.4	3578.5	8	2	C214-213	95	1158.1	2359.3	5	2
							C20-214	99	1538.3	4867.9	11	2

Table 3 :

 3 Cost of changing the cables in the conventional planning

						Year
	Cable	Length,	Old	New	Cost,	to
		[km]	cable	cable	[ke]
						invest
	C2-21	1.15	3x240 3x300 152.9	2
	C213-212	0.17	3x240 3x300 22.3	3
	C214-213	0.12	3x240 3x300 15.6	2
	C20-214	1.77	3x240 3x300 236.5	2

Table 4 :

 4 Cost of the conventional planning solutionAs was previously mentioned, the conventional solution for handling CEs is to invest in the reinforcement of the cables. The costs of changing the cables in the conventional planning and CAPEX, OPEX, TOTEX and VoF are given in the Tables

		Year 3 Year 4
	C conv,CAP EX,y,95 , [ke]	10.1	0.6
	C conv,OP EX,y,95 , [ke]	4.6	0.2
	C conv,T OT EX,y,95 , [ke]	14.7	0.8
	V oF y,95 , [ke]	14.7	0.8/15.5
	6.1. Conventional planning

Table 5 :

 5 Parameters used in the cost functions

	Parameter Value
	Cost of ENS, [e/kVAh]	40
	Cost of losses, [e/kWh]	0.1
	Discount or interest rate, [-]	0.05
	Base cost of energy provided by DR, [e/kWh]	0.1
	Cost of ICT for DR, BESS, [ke/unit]	0.5
	Base cost of energy provided by BESS, [e/kWh]	0.5
	Cost of CB, [ke]	8
	Max. number of SWs, [-] 50000
	Cost of sensors for DLR, [ke/km]	1.0
	Cable LT, [year]	40
	used in Eq. 1 for year 4 is just the cost of the reinforcement of cable C213-212.

Table 6 .

 6 S f lex,y,P % and E f lex,y,P % are used to obtain the costs of using AEs. Due to the network's topology, to avoid any potential CE, an AE should be able to provide S f lex,y,P % and E f lex,y,P % equal to the ones for the main cables C2-21 or C20-214. If the total amount of FSs from any particular type of AEs is not sufficient to cover either S f lex,y,P % or E f lex,y,P % shown in Table2, the maximum capacity and energy that this AE could supply are used for the calculations. The costs are calculated in the assumption that only one type of AE is used in both year 3 and 4. The costs of using AEs in year 4 will depend on whether the DSO decided to utilize their FSs in year 3. If the DSO decided to use AEs in year 3, most of the CAPEX costs would be covered that year. That makes the CAPEX costs in year 4 much smaller as there are fewer new AEs that should be enabled by the DSO. This allows having zero CAPEX for RE in year 4, after it was constructed in year 3.With the exception of RE in year 3, no AE/FS can provide the required capacity and energy alone. Therefore a combination of multiple AEs should be used. The final solution, which uses the combination of DLR and BESS is given in Table7. Comparison of the TOTEX for both years shown in Table7with the VoF in Table4reveals the amount of potential benefits BF f lex for the DSO that could be achieved. Utilizing FSs from DLR and BESS in

	0	100	200	300	400	500	600	700	800	900 1000
				Number of activations,[-]

Table 6 :

 6 Maximum capacity and energy, and costs of each AE/FS available in the distri-

	bution network
		Parameter	DR	BESS	RE	DLR
		S f lex,3,95 , [kVA]	227.5	808.0 1334.4 1284.2
		E f lex,3,95 , [kVAh]	835.3 3288.4 3578.5 3538.5
	Year 3	C f lex,CAP EX,3,95 , [ke]	2.3	1.6	194.3	3.0
		C f lex,OP EX,3,95 , [ke]	1.0	1.5	0.002	0.008
		C f lex,T OT EX,3,95 , [ke]	3.3	3.1	194.3	3.0
		S f lex,4,95 , [kVA]	214.7	808.0 1538.3 1390.6
		E f lex,4,95 , [kVAh]	1173.7 4367.8 4867.9 4742.0
	Year 4	C f lex,CAP EX,4,95 , [ke]	0.9	0.9	0	0.2
		C f lex,OP EX,4,95 , [ke]	1.7	1.4	0.003	0.01
		C f lex,T OT EX,4,95 , [ke]	2.6	2.3	0.003	0.21

Table 7 :

 7 Cost of the solution using flexibility

		Parameter
	Year 3	C f lex,CAP EX,3,95 , [ke]	4.2
	(DLR & BESS )	C f lex,OP EX,3,95 , [ke]	0.1
		C f lex,T OT EX,3,95 , [ke]	4.2
		BF f lex,3,95 , [ke]	10.5
	Year 4	C f lex,CAP EX,4,95 , [ke]	0.2
	(DLR & BESS)	C f lex,OP EX,4,95 , [ke]	0.06
		C f lex,T OT EX,4,95 , [ke]	0.3
		BF f lex,4,95 , [ke]	0.5/15.2

Acknowledgment

This work was supported by the Danish Energy Development Programme (EUDP) through the EnergyLab Nordhavn (Grant : EUDP 64015-0055).

Appendix .1. Initial data

Appendix .3. Modelling of cost functions -DLR

In a traditional operation of distribution networks, the power ratings of overhead lines and cables are considered as fixed values independent of the external conditions. The degree of overloading defined by the coefficient L is typically also set as a constant either for the whole operation period or for a specific season.

However, the exact MCCs of the lines are constantly changing and determined by the ability to dissipate the heat created by an electric current.

Dynamically changing MCC based on the external conditions is referred to as DLR. For the cables buried in the soil, DLR is a function of the soil properties, soil temperature, precipitation and cable's burial depth.

Knowing moisture content, the soil thermal diffusivity δ T can be found:

Soil thermal resistivity:

Eq. A4 shows how the distribution of the soil temperature can be calculated:

The cable temperature can then be found as follows:

Using the cable's actual temperature it is possible to estimate DLR using the following expression from [START_REF] Michiorri | Investigation into the influence of environmental conditions on power system ratings[END_REF]:

Detailed explanations of how to calculate all the parameters shown above could be found in [START_REF] Michiorri | Investigation into the influence of environmental conditions on power system ratings[END_REF][START_REF] Bracale | A new procedure to forecast the dynamic rating of buried cables[END_REF][START_REF] Van Genuchten | A closed-form equation for predicting the hydraulic conductivity of unsaturated soils 1[END_REF].

Since in most of the cases using DLR assumes applying higher than rated electrical stresses to the cable, its LT will be reduced because of that. Such reduction can be used to estimate the cost of DLR via cost function

showing the cost of LT's reduction depending on the duration and power required from DLR during CE:

Using Arrhenius model from Eq. A8 [START_REF] Buhari | Modelling of ageing distribution cable for replacement planning[END_REF] the change in the LT of the cable due to the temperature T cable during ith CE is estimated, with T cable is proportional to I cable and S f lex :

The temperature ratio T is then calculated as:

LT after ith CE referred to the rated temperature LT cable,D f lex,i,end is calculated using Eq. A8 and replacing LT cable,rated with ∆LT cable,i . Finally, ∆LT cable,i :

Using ∆LT cable,i and C LT cable cost function F O1 for DLR is constructed as shown in Fig. A2.


Torbjørn Cunis

Anders La Cour-Harbo

Piecewise Polynomial Model of the Aerodynamic Coefficients of the Cumulus One Unmanned Aircraft

This document illustrates a piecewise polynomial model for Cumulus One, based on continuous fluid dynamics simulation, derived using the pwpfit toolbox. The model has been developed in a joint project of Sky-Watch and the University of Aalborg.

Preliminaries

If not stated otherwise, all variables are in SI units. We will refer to the following axis systems of ISO 1151-1: the body axis system (x f , y f , z f ) aligned with the aircraft's fuselage; the air-path axis system (x a , y a , z a ) defined by the velocity vector V A ; and the normal earth-fixed axis system (x g , y g , z g ). The orientation of the body axes with respect to the normal earth-fixed system is given by the attitude angles Φ, Θ, Ψ and to the air-path system by angle of attack α and side-slip β; the orientation of the air-path axes to the normal earth-fixed system is given by azimuth χ A , inclination γ A , and bank-angle µ A . (Fig. ??.)

x g x a x f z g z a z f mg -Z A F -X A V A Θ γ A α (a) Longitudinal axes (β = µ A = 0). x g x f x a y g y f y a F -X A Y A V A Ψ β χ A
(b) Horizontal axes (γ A = 0).

Fig. 1: Axis systems with angles and vectors. Projections into the plane are marked by .

I. Aerodynamic Coefficients

The piece-wise polynomial models of the aerodynamic coefficients are given

C (α, β, ξ, η, ζ) = C pre (α, β, . . . ) if α ≤ α 0 , C post (α, β, . . . ) else, (1)
where

C ∈ {C X , C Y , C Z , C l , C m , C n }
are polynomials in angle of attack, side-slip, and surface deflections; and the boundary is found at

α 0 = 17.949°. ( 2
)
The polynomials in low and high angle of attack, C pre , C post , are sums

C pre = C pre α (α) + C ξ (β, ξ) + C pre η (α, η) + C ζ (β, ζ) ; (3)
C post = C post α (α) + C ξ (β, ξ) + C post η (α, η) + C ζ (β, ζ) . (4)
In the following subsections, we present the polynomial terms obtained using the pwpfit toolbox. Coefficients of absolute value lower than 10 -2 have been omitted for readability.

A. Domain of low angle of attack

Polynomials in angle of attack:

C pre Xα = -2.566 × 10 -2 + 5.722 × 10 -1 α + 1.496α 2 -1.148 × 10 1 α 3 ;

(5)

C pre Yα = 5.402 × 10 -2 -2.345 × 10 -1 α -2.001α 2 + 7.054α 3 ; (6)
C pre Zα = -3.475 × 10 -1 -5.467α + 1.853α 2 + 2.663 × 10 1 α 3 ; (7)

C pre lα = 4.875 × 10 -2 -2.190 × 10 -1 α -2.004α 2 + 7.146α 3 ; (8)
C pre mα = 6.214 × 10 -2 -1.755α -3.427α 2 + 1.256 × 10 1 α 3 ; (9)

C pre nα = 4.748 × 10 -2 -2.097 × 10 -1 α -2.016α 2 + 7.171α 3 ; (10)
Polynomials in angle of attack and elevator deflections:

C pre Xη = 4.327 × 10 -2 η -4.458 × 10 -1 αη + 3.370 × 10 -1 α 2 η -4.567 × 10 -1 αη 2 + 7.331 × 10 -2 η 3 ; (11)
C pre Yη = 1.832 × 10 -2 η + 7.484 × 10 -2 αη -4.384 × 10 -1 α 2 η; (12)
C pre Zη = -2.567 × 10 -1 η + 3.085 × 10 -1 αη -5.105 × 10 -2 η 2 -7.394 × 10 -1 α 2 η + 6.936 × 10 -1 αη 2 + 1.337 × 10 -1 η 3 ; (13)

C pre lη = 1.699 × 10 -2 η + 8.936 × -2 αη -4.564 × 10 -1 α η + 1.571 × 10 -2 αη 2 ; ( 14
)
C pre mη = -9.028 × 10 -1 η + 7.437 × 10 -1 αη -4.924 × 10 -2 η 2 -8.415 × 10 -1 α 2 η + 2.210αη 2 + 5.251 × 10 -1 η 3 ; (15) C pre nη = 1.532 × 10 -2 η + 8.758 × 10 -2 αη -4.513 × 10 -1 α 2 η + 1.487 × 10 -2 αη 2 . ( 16
)

B. Domain of high angle of attack

Polynomials in angle of attack:

C post Xα = 1.266 × 10 -2 -3.159 × 10 -1 α + 3.832 × 10 -1 α 2 -1.226 × 10 -1 α 3 ; (17)
C post Yα = -2.297 × 10 -2 α 2 + 1.337 × 10 -2 α 3 ; (18)
C post Zα = -4.179 × 10 -1 -2.345α + 9.586 × 10 -1 α 2 -3.665 × 10 -2 α 3 ; (19)
C post lα = 2.006 × 10 -2 -8.200 × 10 -2 α + 1.012 × 10 -1 α 2 -3.515 × 10 -2 α 3 ; (20) C post mα = -2.552 × 10 -1 -5.131 × 10 -1 α -2.677 × 10 -1 α 2 + 1.332 × 10 -1 α 3 ; (21) C post nα = 1.159 × 10 -2 -3.062 × 10 -2 α + 2.727 × 10 -2 α 2 ; (22)
Polynomials in angle of attack and elevator deflections:

C post Xη = -1.342 × 10 -2 η -1.663 × 10 -1 αη -1.796 × 10 -1 η 2 + 2.254 × 10 -2 α 2 η + 9.274 × 10 -2 αη 2 + 7.331 × 10 -2 η 3 ; (23) C post Zη = -2.922 × 10 -1 η + 1.829 × 10 -1 αη + 1.277 × 10 -1 η 2 + 2.284 × 10 -2 α 2 η + 1.229 × 10 -1 αη 2 + 1.337 × 10 -1 η 3 ; (24) C post mη = -9.498 × 10 -1 η + 6.099 × 10 -1 αη + 5.093 × 10 -1 η 2 + 6.456 × 10 -2 α 2 η + 4.264 × 10 -1 αη 2 + 5.251 × 10 -1 η 3 ; (25) C post Yη , C post lη , C post nη are approximately zero. ( 26
)
C. Full-envelope polynomials Polynomials in side-slip and aileron deflections:

C Xξ = 6.557 × 10 -2 β 2 + 4.214 × 10 -2 βξ -1.493ξ 2 + 4.264 × 10 -2 ξ 3 -1.647 × 10 -2 β 4 -2.321 × 10 -2 β 3 ξ + 9.649 × 10 -2 β 2 ξ 2 + 2.859 × 10 -2 βξ 3 + 3.084 × 10 1 ξ 4 ;
(27)

C Yξ = -3.697 × 10 -1 β -1.570 × 10 -1 ξ -3.231 × 10 -2 β 2 -6.137ξ 2 + 7.416 × 10 -2 β 3 + 1.611 × 10 -2 βξ 2 + 2.214ξ 3 + 6.487 × 10 -1 β 2 ξ 2 + 5.323 × 10 -2 βξ 3 + 1.456 × 10 2 ξ 4 ;
(28)

C Zξ =
3.411 × 10 -1 β 2 + 4.141 × 10 -1 βξ -3.717ξ 2 + 4.234 × 10 -2 ξ 3 -9.915 × 10 -2 β -1.922 × 10 -1 β 3 ξ + 1.945 × 10 -1 β 2 ξ 2 + 9.909 × 10 -1 βξ 3 + 9.856 × 10 1 ξ 4 ;

(29)

C lξ = -5.798 × 10 -2 β -3.929 × 10 -1 ξ -2.906 × 10 -2 β 2 -5.551ξ 2 + 1.763 × 10 -2 β 3 + 1.722 × 10 -1 β 2 ξ + 4.557 × 10 -1 ξ 3 + 5.835 × 10 -1 β 2 ξ 2 + 5.332 × 10 -2 βξ 3 + 1.319 × 10 2 ξ 4 ;
(30)

C mξ = -3.978 × 10 -2 β 2 + 5.554 × 10 -1 βξ -4.689ξ 2 + 4.229 × 10 -2 ξ 3 + 2.585 × 10 -2 β 4 -2.309 × 10 -1 β 3 ξ + 4.077 × 10 -1 β 2 ξ 2 -5.068 × 10 -1 βξ 3 + 1.160 × 10 2 ξ 4 ; (31) C nξ = 3.686 × 10 -2 β + 1.392 × 10 -2 ξ -2.828 × 10 -2 β 2 -5.410ξ 2 -1.160 × 10 -1 ξ 3 + 5.678 × 10 -1 β 2 ξ 2 + 5.334 × 10 -2 βξ 3 + 1.286 × 10 2 ξ 4 . ( 32
)
Polynomials in side-slip and rudder deflections:

C Xζ = -1.031 × 10 -2 βζ -7.986 × 10 -2 ζ 2 + 1.086 × 10 -2 βζ 3 + 9.739 × 10 -2 ζ 4 ; (33)
C Yζ = -8.526 × 10 -2 ζ -2.128 × 10 -2 β 2 -3.569 × 10 -1 ζ 2 + 1.049 × 10 -2 β 2 ζ + 1.924 × 10 -2 βζ 2 + 7.679 × 10 -2 ζ 3 + 5.663 × 10 -2 β 2 ζ 2 + 5.247 × 10 -1 ζ 4 ;
(34)

C Zζ = -2.036 × 10 -2 β 2 -1.332 × 10 -2 βζ -1.743 × 10 -1 ζ 2 + 1.834 × 10 -2 β 2 ζ 2 + 3.683 × 10 -2 βζ 3 + 2.631 × 10 -1 ζ 4 ; (35)
C lζ = -1.930 × 10 -2 β 2 -3.221 × 10 -1 ζ 2 + 5.111 × 10 -2 β 2 ζ 2 + 4.735 × 10 -1 ζ 4 ; (36)
C mζ = -3.231 × 10 -2 β 2 -1.436 × 10 -2 βζ -2.672 × 10 -1 ζ 2 + 1.241 × 10 -2 β 3 ζ + 1.472 × 10 -2 β 2 ζ 2 + 1.051 × 10 -2 βζ 3 + 4.193 × 10 -1 ζ 4 ;
(37)

C nζ = 2.158 × 10 -2 ζ -1.882 × 10 -2 β 2 -3.137 × 10 -1 ζ 2 -1.338 × 10 -2 ζ 3 + 4.978 × 10 -2 β 2 ζ 2 + 4.611 × 10 -1 ζ 4 . ( 38
)

II. Equations of Motion

We have the aerodynamic forces and moments in body axis system

  X A Y A Z A   f = 1 2 SV 2 A   C X C Y C Z   ;   L A M A N A   f = 1 2 SV 2 A   b C l c A C m b C n   +   X A Y A Z A   f × x cg -x ref cg ; (39)
and the weight force by rotation into the body axis system

X G g = -g sin Θ;
(40)

Y G g = -g sin Φ cos Θ; (41) Z G g = -g cos Φ cos Θ; (42)
The resulting forces lead to changes in the velocity vector, in body axis system, by

VAf = 1 m   X A + X G + X F Y A + Y G Z A + Z G   f +   p q r   × V Af . ( 43
)
with the thrust X F f = F (engines aligned with the x f -axis). For a symmetric plane (I xy = I yz = 0), the resulting moments in body axis are given as

L f = L A f -q r (I z -I y ) + p q I zx ; (44) M f = M A f -p r (I x -I z ) -p 2 -r 2 I zx ; (45) N f = N A f -p q (I y -I x ) -q r I zx ; (46)
with the inertias I x , I y , I z , I zx . The changes of angular body rates are then given as

ṗ = 1 I x I z -I 2 zx (I z L f + I zx N f ) ; (47) q = 1 I y M f ; (48) ṙ = 1 I x I z -I 2 zx zx L f + I x N f ) . (49)
Here, the normalized body rates have been used with

  p q r  = 1 2V A   b p c A q b r   ⇐⇒   p q r   = 2V A   b -1 p c -1 A q b -1 r  (50)
and

  ṗ q ṙ  = 1 V A   1 2   b ṗ c A q b ṙ   -VA   p q r    . (51)
The change of attitude is finally obtained by rotation into normal earth-fixed axis system:

Φ = p + q sin Φ tan Θ + r cos Φ tan Θ; (52) Θ = q cos Φ -r sin Φ;
(53) Ψ = q sin Φ cos -1 Θ + r cos Φ cos -1 Θ.

(54)

III. Longitudinal Model

The longitudinal model is restricted to the x a -z a -plane, assuming β = µ A = χ A = 0.

A. Aerodynamic Coefficients

The longitudinal aerodynamic coefficients, C L , C D , C m , are obtained as

C L C D = + sin α -cos α -cos α -sin α C X C Z . (55)

B. Equations of motion

The longitudinal equations of motion are given as x a , y a , z a = Air-path axis system; x f , y f , z f = Body axis system; x g , y g , z g = Normal earth-fixed axis system;

VA = 1 m T cos α - 1 2 ρSV 2 A C D (α, η, q) -mg sin γ A , (56)
γA = 1 mV A T sin α + 1 2 ρSV 2 A C L (α, η, q) -mg cos γ A , (57)
q = 1 I y 1 2 ρSc A V 2 A C m (α, η, q) ; (58

  angle of attack boundary (°); β = Side-slip angle (rad); γ A = Flight-path angle relative to air (rad); ζ = Rudder deflection (rad), negative if leading to positive yaw moment; η = Elevator deflection (rad), negative if leading to positive pitch moment; Θ = Pitch angle (rad); µ A = Air-path bank angle (rad); ξ = Aileron deflection (rad), negative if leading to positive roll moment; = Air density ( = 1.200 kg m -3 ); χ A = Air-path azimuth angle (rad); Φ = Bank angle (rad); Ψ = Azimuth angle (rad); C l = Aerodynamic coefficient moment body x f -axis (•); C m = Aerodynamic coefficient moment body y f -axis (•); C n = Aerodynamic coefficient moment body z f -axis (•); C D = Aerodynamic drag coefficient (•), positive along negative air-path x a -axis; C L = Aerodynamic lift coefficient (•), positive along negative air-path z a -axis; C X = Aerodynamic coefficient force body x f -axis (•); C Y = Aerodynamic coefficient force body y f -axis (•); C Z = Aerodynamic coefficient force body z f -axis (•); D = Drag force (N), positive along negative air-path x a -axis; F = Thrust force (N), positive along body x f -axis; L = Lift force (N), positive along negative air-path z a -axis; q = Pitch rate (rad s -1 ); b = Reference aerodynamic span (b = 2.088 m); c A = Aerodynamic mean chord (c A = 2.800 × 10 -1 m); g = Standard gravitational acceleration (g ≈ 9.810 m s -2 ); m = Aircraft mass (m = 2.619 × 10 1 kg); p = Roll rate (rad s -1 ); r = Yaw rate (rad s -1 ); L f = Roll moment (N m), mathematically positive around x f -axis; M f = Pitch moment (N m), mathematically positive around y f -axis; N f = Yaw moment (N m), mathematically positive around z f -axis; S = Wing area (S = 5.500 × 10 -1 m 2 ); V A , V A = Aircraft speed and velocity relative to air (V A = V A 2 , m s -1 ); X A f = Aerodynamic force along body x f -axis (N); X F f = Thrust force along body x f -axis (N); Y A f = Aerodynamic force along body y f -axis (N); Z A f = Aerodynamic force along body z f -axis (N); (•) post = Domain of high angle of attack; (•) pre = Domain of low angle of attack;

ONERA -The French Aerospace Lab, Department of Information Processing and Systems, Centre Midi-Pyrénées, Toulouse, 31055, France, e-mail: torbjoern.cunis@onera.fr; and ENAC, Université de Toulouse, Drones Research Group, Toulouse, 31055, France.

Aalborg University, Department of Electronic Systems, Aalborg East, 9220, Denmark, e-mail: alc@es.aau.dk.


Rémi Guidoux

Gael Manautines

Arnaud Dion

Jean-Baptiste Munier

Vincent Calmettes
email: vincent.calmettes@supaero.fr

Real Time Estimation of rigid body orientation based on inertial sensors measurements

This paper presents a programmable DSP board developed for real-time estimation of rigid body orientation. This work takes place in a research project carried out by students. The project consists of unmanned aerial vehicle's (UAV) attitude estimation. This estimation is processed by an Inertial Reference System (IRS) including the TMS320F2812 processor, low cost integrated Micro Electro Mechanical Systems (iMEMS) such as accelerometers and gyroscopes, temperature sensors, and magnetometers. Gyroscopes provide angular velocities that can be integrated to yield orientation. Accelerometers and magnetometers, pointing to gravity direction and north direction enables drift compensation. So the system performs a drift-free attitude estimation using quaternions for rotations representation. One benefit of this system is that it enables high integration by exploiting the processor's capacities for sensors interfacing. Moreover it takes advantage of the high performance processor to achieve a real time estimation of the body's orientation. For that digital processing includes signals filtering, measurements compensation and a lot of trigonometric functions for Euler angles computing. The system performances can be measured using a PC for results visualisation. The system has been designed to integrate in the future GPS component and pressure sensor. These components will allow more accurate attitude estimation and positioning.

Introduction

Supaero has been working on unmanned aerial vehicle (UAV) development for several years. UAV can be piloted by high level orders provided by guidance unit and measurements provided by navigation system. Navigation system is used for attitude and position estimation. Attitude estimation consists of accurate measurements of yaw, roll and pitch that are required for determining a reference frame. Positioning is the expression of the vehicle location in this reference frame. Such a system usually uses GPS coupled to an Inertial Reference System (IRS) [START_REF] Shin | A Quaternion-based unscented Kalman filter for the integration of GPS and MEMS INS[END_REF]. IRS provides measurements with a good short term accuracy and doesn't depend on external signal. However these measurements are affected by low drift. So this system must be calibrated, using GPS measurements for example. Nevertheless, in certain environment such as urban environment, GPS measurements can be corrupted at any moments during long time period. When a low cost IRS is used auxiliary measurements are necessary to provide drift-free attitude. Especially accelerometers and magnetometers can be used to define a North, Down, East frame [START_REF] Luis Marins | An extended Kalman filter for Quaternion-Based Orientation Estimation Using MARG Sensors[END_REF].

The aim of this project is to develop a compact drift-free IRS for attitude estimation. This low-cost onboard unit has been developed by students. It includes low cost sensors that provide measurements processed by the TMS 320F2812 DSP. System integration is described, particularly the major components and DSP. Then the estimation methodology is presented by insisting on the main operations that affect the computing load. Finally the final product is characterised and extension is proposed allowing GPS/IRS hybridisation and so vehicle positioning.

System architecture

This IRS system consists of a compact unit based on the TMS320F2812 processor. This unit is designed to yield 3D measurements by integrating sensors mounted on 3 orthogonal boards. Tri-axial Gyroscope provides angular velocities around three orthogonal axis, which must be integrated to estimate orientation. To compensate orientation drift, tri-axial accelerometer and tri-axial magnetometers are implanted. These components permit to determine respectively the gravity and the north direction. Moreover measurements provided by the accelerometers could be used in the final project by the navigation algorithm. All the measurements are temperature sensitive. Compensation needs to measure the temperature with a temperature sensor. Finally the system integrates a barometric pressure sensor and a GPS component. These components are not used by the algorithm presented in this paper. In the final project they will enable attitude estimation improvement and position estimation.

Gyroscopes

Gyroscope is used to achieve angular motion measurement, providing an output signal proportional to the angular rate. A monolithic integrated circuit from Analog Devices (ADXRS300) is proposed for this system. This low cost component combines mechanical structures with signal processing circuitry. Three identical components are used for measurements around three orthogonal axis. Signals delivered by the gyroscopes are amplified and applied to the processor's ADC module. The main features of the gyroscope are given here:

-Range : ±300 °/s -Sensitivity : 5mV/(°/s) -Bandwidth : DC to 2 KHz -Noise density : 0.1 (°/s)/√Hz

Accelerometers

By considering that the acceleration is small and effective for a short duration, accelerometers can be used to determine gravity direction. For this application iMEMS accelerometers from Analog Devices (ADXL202JE) are implanted. This component includes 2 accelerometers and signal processing circuitry. It provides 2 PWM signals at 1KHz, with duty cycles depending on the acceleration along 2 orthogonal axis. PWM conditioning enables direct connection to the capture unit of the processor's event manager. So this event manager is used for duty cycle measurement. The main features of the accelerometer are given here:

-Range : [-2g ; +2g] -Sensitivity : 12.5 %/g (PWM Duty Cycle) -Bandwidth : 6 KHz -Noise density : 200µg√Hz

Magnetometers

Magnetic sensor achieves measurement of the magnitude of external magnetic fields. So three sensors, by sensing the Earth's magnetic field, can be used as heading sensor. For this system a solid state component has been selected. This component delivered by Honeywell provides a differential signal proportional to the magnetic field applied along its main axis. Signal conditioning such as Differential to common mode conversion is added at the input of the TMS320F2812 ADC module. The main features of the magnetic sensor are given here :

-Sensitivity depending on the power level : 1mV/V/Gauss -Noise density: 16. 10 -6 Gauss√Hz -Bandwidth : 5MHz -Heading error taking into account noise, hysteresis and linearity : 0.15°

Temperature sensor

The Maxim MAX6666 temperature sensor is used for IRS sensors calibration (temperature compensation). Providing a PWM signal at 35Hz it can be directly connected to the capture unit. As it is done for accelerometers, the processor's event manager is used for duty cycle measurement.

Global architecture

The global architecture is given on the Figure 1. [ ]

) ( . ref T T T k b M K P - - - = (1)
P is the corrected magnitude vector when M is the measured magnitude vector, b vector represents bias measurements, T is the measured temperature and T ref the reference temperature, k T is a scale factor depending on the sensitivity to the temperature, and K is a matrix deduced from the sensor sensitivity and from the geometry.

Attitude Estimation

Attitude can be defined by a set of angular positions from a frame to another reference frame.

Here the earth-fixed reference frame is the North-East-Down coordinates frame centred at the initial position of the body (NED frame). To express the rotation from the body frame to this reference frame, different techniques can be used (Euler angles, quaternions, rotation matrix).

For this application quaternions are preferred to Euler angles to eliminate singularities associated with attitude estimation [START_REF] Shin | A Quaternion-based unscented Kalman filter for the integration of GPS and MEMS INS[END_REF], [START_REF] Luis Marins | An extended Kalman filter for Quaternion-Based Orientation Estimation Using MARG Sensors[END_REF]. So equations of rotational motion are integrated using quaternions integrator. The rotational matrix is deduced from the angular rates provided by gyroscopes. Due to the low accuracy of such sensors (gyro drift of about 1000 deg/h), the attitude estimated from these measurements suffers errors. Therefore attitude correction is applied to achieve improved performance. Figure 3 represents the software actigram.

Initialisation

The system must be initialised to determine the initial values of the quaternions. Measurements provided by tri-axial accelerometer and tri-axial magnetometer are averaged during 3 sec and used to determine the Gravity direction and the North direction in the mobile frame, and so to express the coordinates of a North, East, Down frame ( )

m m m D E N , ,
in the mobile frame. These coordinates are used for Euler angles computing : And initial quaternions init Q can be expressed from these angles :

      Ψ             -       Ψ             =       Ψ             +       Ψ             =       Ψ             -       Ψ             =       Ψ             +       Ψ            θ φ θ φ q q q q (3)

Quaternions Integration

Tri-axial gyroscope provide angular rates (p, q, r) for rotation estimation. This rotation can be expressed by updating the quaternions vector.

Φ, θ, Ψ c Q Q init Q dQ
Quaternions integration is performed using the Wilcox method to integrate the differential equation :

= = + =             =             - - - - - - - = = + • e e e T T A Q A X X T I X Q q q q q p q r p r q q r p r q p AQ Q (4)

Quaternions Bias Estimation

Quaternions, which then define the roll, pitch and yaw angles, are used to compute the transformation matrix from the reference frame(NED frame) to the mobile frame. So the theoretical Gravity vector, oriented towards the Earth centre can be expressed in the mobile frame, using this transformation matrix. It can be compared to local vertical coordinates returned by the tri-axial accelerometer. In the same way, local North vector is defined in the mobile frame and compared to the local magnetic field provided by the tri-axial magnetometer. The differences on down and north coordinates are induced by errors on quaternions estimates. These errors can be computed by taking advantage of the orthogonal projection properties of the pseudo inverse.

The transformation matrix from the reference frame to the mobile frame is defined from quaternions

c Q : ( ) ( ) ( ) ( ) ( ) ( )           + - - - + + - + - - + + - - + = 2 3 2 2 2 1 2 0 1 0 3 2 2 0 3 1 1 0 3 2 2 3 2 2 2 1 2 0 3 0 2 1 2 0 3 1 3 0 2 1 2 3 2 2 2 1 2 0 2 2 2 2 2
2 2 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q R m N (5)

The Gravity vector expressed in the reference frame is :

T NED g D ] 0 0 [ = (6)
North magnitude (m) and inclination (I) depend on the location. These parameters have to be determined to express the North vector :

( ) ( ) T T NED I I m n n m N ] cos 0 [sin . ] 0 [ . 2 1 = = (7)
Using the transformation matrix these vectors can be defined in the mobile frame.

NED m N m NED m N m N R N D R D 2 2 , = = (8)
An error vector, resulting from bias on quaternions estimate, is obtained by comparing the coordinates of these vectors to those of the vectors B and A which represent tri-axial magnetometer and accelerometer measurements.

A D E B N E E E E m G m N N G - = - =         = and where , (9)
By noting the Jacobian Matrix X (6 raws, 4 columns), biases on quaternions dQ are expressed from this error vector E :

( )

            ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ = = - 3 2 1 0 3 2 1 0 1 X where * q m N q m N q m N q m N q m D q m D q m D q m D E X X X dQ T T (10)
Greville Algorithm [START_REF] Greville | The pseudo inverse of a rectangular or singular matrix and its applications to the solution of system of linear equations[END_REF] has been implemented to compute the pseudo-inverse required by this processing.

Quaternions Correction

The error vector dQ is used to correct quaternions after integration. The corrector model is given here :

Euler angles computing

Orientation is defined by three angles : Pitch and Roll. These angles can be deduced from the quaternions :

( ) ( ) ( ) ( ) ( ) ( ) 2 3 2 2 2 1 2 0 3 2 1 0 3 1 2 0 2 3 2 2 2 1 2 0 2 1 3 0 , 2 2 tan , 2 sin , , 2 2
tan q q q q q q q q a q q q q a q q q q q q q q a

+ - - + = - = - - + + = Ψ φ θ (11)
After computing Yaw, Pitch and Roll are transmitted by SPI interface to a PC. Used during the development step, PC allows data visualisation and IRS assessment.

Assessment

Figure 5 shows the system structure. The code is written in C. Special configuration enables hardware validation. In this mode, PC connection displays samples of filtered signals. In other respects, the filtering algorithm was executed on DSP, using pre-registered data. It achieves good estimation. The algorithm uses a lot of trigonometric functions. Data have been represented in IQ24 format and the IQmath library permits to produce highly efficient assembly code. Finally the code has been applied to sensors samples for yielding attitude estimation. In that case the system provides Yaw, Pitch and Roll. These angles are transferred to the PC for visualisation. Up to now the system couldn't be calibrated. Especially

Integration Wilcox Method z -1 k P k I + + + + + - c Q Q dQ
Gyroscopes pre-processing requires calibration to determine magnitude and orthogonality corrections. Appropriate test bench has to be used to apply known angular rates.

The length of the Flash memory used for code implementation is 1C0A words. The computing load has been evaluated by measuring the duration of the processing. This duration is 200uS or 50uS depending on whether the code is executed in flash or internal Ram memory. These values are largely lower than the sample period (T e =10mS) and this application doesn't need to copy time-critical portions of code from Flash memory into internal RAM memory.

Conclusion and future prospect

Students have contributed to this project. Filtering methodology, including Kalman filtering for GPS hybridization has been proposed and studied by a group of 4 second-year students from SUPAERO, during an internal project of one month. A second-year student, from an Electronic Engineering School (ESME) has carried on this study during a summer project, with the aim to develop the C code for the IRS module. The present board was designed by a student from IPST (Engineering school for adults) during his final project. Algorithms concerning the IRS were implanted by another student from IPST. The workload, including study, coding, hardware conception, represents approximately one student year. The next step of this project concerns calibration. Finally GPS/IRS hybridization is to considered. This implementation requires algorithm coding and evaluation in term of computing load.

Figure 1 :

 1 Figure 1: IRS architecture

Figure 2 :

 2 Figure 2: Signals filtering. 10 phases-FIR Filter

Figure 3 :

 3 Figure 3: Software Actigram

Figure 4 :

 4 Figure 4: Control unit

Figure 5 :

 5 Figure 5: System integration

  Heading signals and angular velocities, provided by magnetic sensors and gyroscopes, are represented in a 12 bits format at the output of the ADC module. Accelerations and temperature are determined by measuring the duty cycle of the PWM signals with the event managers. Event managers registers permit a 16 bits representation.

	3 Pre-processing
	Sensors signals (accelerations, heading signals, angular velocities) are sampled at 1KHz but
	processed for attitude estimation at the rate of 100 Hz. A FIR filter is used for sensors signals
	denoising. Filtering is followed by decimation for sample rate conversion. These two
	operations can be processed simultaneously using a 50 coefficients-polyphase FIR filter.
	Including 10 branches, this architecture reduces the computing load to 5 operations
	(multiplication-accumulation) per mS and per channel.
		1ms	z -1	Delays: 10mS	z -1
			h 0		h 10		h 40
	Sensor signal
	Fe 1 =1KHz	h 1 z -1		h 11	z -1	h 41	+	Filtered signal Fe 2 =1KHz
			z -1			z -1
			h 9		h 19		h 49
		X, Y	Capture
		accelerometers	unit
		Temperature
		Sensor
		X, Z accelerometers	Capture unit	EV B	Serial Interface	UART	Test
				Timer 4	SCI
		Differential
	X, Y, Z	to
	magnetometers	common
		Conversions
	X, Y, Z velocity rate	Amplifiers	ADC module	Serial Interface SCI	GPS Component
	Pressure
	Sensor