Abdelouadoud Seddik Yassine
email: seddik.abdelouadoud@mines-paristech.fr

Robin Girard

François-Pascal Neirac

Guiot Thierry
email: thierry.guiot@cstb.fr

A

Keywords: Distributed power generation, power system planning, smart grids I. CONTEXT

As the interest for distributed generation and storage grows, so does the need for new tools to assist in the planning of the distribution network. Indeed, as the current passive network transforms into an Active Distribution Network (ADN) [START_REF] Pilo | Active distribution network evolution in different regulatory environments[END_REF] with the introduction of partially and totally controllable generation and storage means, planning studies based solely on power flows for extreme load conditions will not be adapted anymore. Considering the similarities between the current transmission network and the future ADN, it is a safe bet to assume that the Optimal Power Flow, a tool first introduced in 1962 by Carpentier [START_REF] Carpentier | Contribution à l'étude du dispatching économique[END_REF] and now widely used for the planning and operation of the transmission network, will prove useful for this purpose. Consequently, the adaptation of the OPF concept and resolution algorithms to the distribution network has been the subject of numerous publications in the last decade such as [START_REF] Gabash | Active-Reactive Optimal Power Flow in Distribution Networks With Embedded Generation and Battery Storage[END_REF], [START_REF] Ahmadi | Optimal power flow for autonomous regional active network management system[END_REF], [START_REF] Swarnkar | Optimal power flow of large distribution system solution for combined economic emission dispatch problem using partical swarm optimization[END_REF], [START_REF] Dolan | Using optimal power flow for management of power flows in active distribution networks within thermal constraints[END_REF]. Among the structural differences between transmission and distribution network, the radial topology and higher R/X ratio of line impedances are often identified as the main obstacles to the direct use of transmission network OPF algorithms to distribution systems as explained in [START_REF] Bruno | Unbalanced three-phase optimal power flow for smart grids[END_REF]. With this in mind, an OPF for a distribution system will still be a highly dimensional nonlinear nonconvex programming problem. Thus, in a planning study for which OPFs have to be solved for each hour of the year and multiple planning options evaluated, the computational burden can become dissuasive. And so, heuristics that allow to decrease the computational burden while controlling the loss of optimality can have practical applications, as demonstrated in [START_REF] Paudyaly | Three-phase distribution OPF in smart grids: Optimality versus computational burden[END_REF]. In the work presented here, we make use of the characteristics of the problem studied to propose a methodology aiming at evaluating the voltage constraints only when they are binding, thus drastically reducing the computational cost of the resolution.

II. METHODOLOGY

A. Problem studied

First, let us define the general form of the OPF problem for which the methodology presented here applies. It consists in finding, for each time steps, P, Q, V and α that verify:

) ( min arg P f P (1) 0 ) , , , (   Q P V DF (2) 0 ) , (   V h (3) 0 ) , (  Q P k (4) 0 ) , (  Q P m (5)
where P and Q are the vectors of, respectively, controllable active and reactive power algebraic injections, V is the vector of node voltage magnitude, α is the tap setting of the PowerTech 2013, Towards carbon free society through smarter grids, 16-20 June 2013, Grenoble, France substation transformer on-load tap changer (OLTC), f represents the objective function, DF the set of DistFlow network equations used in, for example, [START_REF] Dukpa | An accurate voltage solution method of radial distribution system[END_REF], h the inequality voltage constraints and k and m represent the other equality and inequality constraints implemented in the problem (e.g. limit on apparent powers, reserve provision requirements or limits on storage means energy capacity). We emphasize here the fact that the methodology proposed is applicable only if the reactive power injections have no influence on the objective function.

B. Problem simplification

In many instances, the problem complexity can be greatly reduced by ignoring ( 2) and (3). Indeed, for example, if the objective function consists in the sum of generation piecewise linear costs and (5) represents bounds on active and reactive power, the problem defined by ( 1), ( 4) and ( 5) is a simple linear programming problem. Moreover, we aim at taking advantage of the fact that voltage limit violation occurrence are infrequent phenomena, as evidenced by [START_REF] Widén | Impacts of distributed photovoltaics on network voltages: Stochastic simulations of three Swedish low-voltage distribution grids[END_REF]. This finding supports the assumption that (2) and (3) will be binding constraints only during a relatively small number of time steps. To support this claim, we present the results of a simulation undertaken on a 69-bus medium-voltage distribution system (see [START_REF] Baran | Optimal capacitor placement on radial distribution system[END_REF] for the detailed characteristics of the network). Residential-type loads are simulated for each bus and scaled so that the maximal annual active load is equal to the one used in [START_REF] Baran | Optimal capacitor placement on radial distribution system[END_REF] and the power factor remains constant. Then, 17 buses are equipped with photovoltaic and energy storage systems in increasing order of driving-point impedance magnitude, with varying penetration levels (the penetration level at a given bus is defined by the ratio between the maximum power of the storage and photovoltaic systems and the annual maximal load at this bus). The storage systems have a 5 hour discharge time and their charge and discharge efficiency are both equal to 0.95. The objective function is defined as a measure of market transaction cost. In the graph below, we observe the number of time steps where the voltage constraints are active or binding, as a function of the penetration level. Consequently, we propose to separate the initial problem into a master problem composed of ( 1), ( 4) and ( 5) that is easier to solve and a slave problem consisting in finding P,Q,V, θ and α that verify:

  2 ) ( min arg P P master P , ( 7
)
while satisfying the constraints ( 2), ( 3), ( 4) and ( 5). The slave problem is still a highly dimensional nonlinear nonconvex programming problem. We will now introduce a procedure that will allow us to solve it only when it is necessary, i.e. when its result has an impact on the objective function of the initial problem.

C. Criticality Criterion 1) Guiding principle

For the remainder of this paper, we will adopt the following convention: a time step will be deemed critical if (3) is a binding constraint, semi-critical if it is active but not binding, and noncritical if it is inactive. On a more practical level, the noncritical time steps are those for which the active powers resulting from the master problem directly verify (3), the semi-critical are those for which the reactive power and/or tap settings but not the active powers have to be adjusted and the critical time steps are those for which the active powers have to be modified. In order to solve the slave problem only when necessary, we need a reliable and computationally effective method to separate the critical time steps from the others. To achieve this, we calculate a relative criticality criterion aimed at ranking the time steps so that the critical time steps are first, followed by the semi-critical and noncritical ones. First, let us set aside the non-critical time steps as these can be identified with a mere power flow to concentrate on separating the semi-critical from the critical ones. Essentially, this means finding a measure of the ability of the reactive powers and OLTC to correct the voltage violations incurred from the behavior of the loads, distributed generation and storage systems resulting from the master problem solution.

2

) Voltage violation and OLTC influence

In order to take into account the influence of the on-load tap changer, we define  as the set of permissible tap-settings and t OLTC V ,  the substation downstream voltage at time step t corresponding to the tap position  . The first step of the criticality criterion calculation is to obtain a power flow for a chosen initial tap setting 0  and then expand the result to all possible tap settings for each bus k by making the following approximation:

.

t OLTC t OLTC t k t k V V V V , , , , 0 0        (8)
We then obtain an estimation of the voltage limit violation during time step t at bus k for the tap setting  in the following way:

             min , , min max , min max , max , , if , if , 0 if , V V V V V V V V V V V V t k t k t k t k t k t k       . ( 9
)
This gives us an estimation of the absolute voltage violation to be compensated at each bus and for each permissible tap-setting.

3) Reactive power voltage violation compensation capability

We base our reasoning on the concept of Voltage Change Potential introduced in [START_REF] Auld | Analysis and visualization method for understanding the voltage effect of distributed energy resources on the electric power system[END_REF] and aimed at evaluating the local voltage variations caused by the connection of a given distributed energy resource. It is based on a linearization of the network equations, under the assumption that the actual voltage change is small compared to the pre-existing voltage, and is defined as follows:

V P R Q X V      , ( 10
)
where X is the driving-point reactance, R the driving-point resistance, Q the reactive power injected at the considered bus, P the active power and V the pre-existing voltage.

For our purpose, we define the Reactive Voltage Change Potential during time step t at bus k as follows:

t k t k k t k Q V Q X V , , ,    (11) t k Q V , 
is effectively a measure of the ability of a given reactive power source to modify the voltage at its bus.

4) Formal definition of the criterion

We now have at our disposal a measure of the voltage violations taking into account the role of the OLTC and an estimation of the reactive power local compensation capability. We can combine them in order to define a local compensated voltage violation in the following way:

t k Q t k t k Q V V V , , , ,        (12)
As we need a single measurement by time step in order to establish our criticality criterion, we define it as the minimum of the L1-norm of the compensated voltage violation, while respecting constraints on reactive power availability and permissible tap-settings. This can be formally expressed as:

    K k t k Q Q t V C , , , min   (13) t k t k Q Q , max , 0   (14)    (15)
Where K is the set of network buses and t k Q , max is the available reactive power at bus k during time step t .

5) Calculation

As  is a discrete variable, we can define, for each    :

    K k t k Q Q t V C , , min   , ( 16
)
while respecting (14). We can observe that, for a fixed  , the terms in the sum are independent from one another. Thus, we also have :

    K k t k Q Q t V C , , min   (17)
where :

          otherwise 0 0 if min , , , , , , max max t k Q t k Q t k Q Q V V V    (18)
We can thus easily calculate

III. APPLICATION

A. Case Study 1) Criticality criterion calculation

We apply this methodology to the case study introduced in paragraph B. The master problem is solved using IBM ILOG CPLEX and the slave problem is solved using IPOPT [START_REF] Wächter | On the Implementation of a Primal-Dual Interior Point Filter Line Search Algorithm for Large-Scale Nonlinear Programming[END_REF], both through their Matlab interface. Below are two graphs presenting the ranking obtained with the criticality criterion for the various types of time steps. In the first simulation considered, 17 buses have been equipped with photovoltaic generators and storage systems, and the level of penetration has been set at a 100%, while in the second, 22 buses have been equipped with a level of penetration of 110%. We can observe that the criticality criterion fulfills its goal, as we obtained a strict separation between the critical time steps and the semi-critical ones.

2

) Computational time spared

As the stated goal of this methodology is to decrease the computational burden associated with such calculations, we now present an evaluation of the computational time spared for various levels of penetration and in the case where 17 buses have been equipped. These results have been on obtained on a desktop computer with a 2.4 GHz processor and 6 GB of RAM. We observe that the computational gain obtained increase substantially starting from a level of penetration of 60%, which also corresponds to the occurrence of a significant and growing number of semi-critical time steps (see fig. 1). This is to be expected, as the slave problem takes typically 0.02 seconds to be solved for non-critical time steps, while it takes between 2.5 and 3 seconds for a semi-critical time step.

3) Voltage control impact

A possible application of the methodology introduced is to evaluate the deoptimization introduced to the master problem scheduling by the voltage control implemented. First, let us observe the solution from the slave problem for a critical time step. We can see that the active power injections from the storage systems are effectively reduced to cope with overvoltage occurrences. This has two effects: the first one is the decrease in voltage due to lower active power injections and the second one is the freeing up of reactive power capabilities that allows for reactive power consumption, further reducing the voltage.

To evaluate the impact of the voltage control, we first compute the gains obtained in terms of scheduling costs through the deployment of photovoltaic generators associated with storage systems, for various scenarios. These are expressed in percentage of the cost of scheduling when only the loads are present. Then we compute the cost associated with voltage control in terms of deoptimization of the scheduling, and in this case, the results are also expressed in percentage of the cost of scheduling when only the loads are present.

B. Limitation of the criticality criterion

We have presented results derived from the use of the criticality criterion in cases where the separation obtained was satisfactory. However, this is not always the case, especially when the level of penetration is increased substantially. For example, the ranking obtained when 22 nodes are equipped and the level of penetration is 140 percents is as follows: The ranking obtained is not satisfactory, as a significant number of critical time steps have lower criticality criterion values than some semi-critical time steps (hereafter named missed critical time steps). This is to be expected, in the sense that the criticality criterion, in essence, aims at describing the behaviors of heterogeneous buses with only one value. Thus, when the behaviors, for a given time steps, of various buses starts to diverge significantly (for example, a high load is present on a bus without storage while a high level of power is injected from another bus), we cannot expect the criticality criterion to be efficient. In order to investigate the limits to the efficiency of the criticality criterion, we define the maximal difference of active voltage change potential, as a proxy measurement of the divergence in behaviors between various buses. It is calculated in the following way: We observe here that the criticality criterion remains efficient until a level of penetration of 130 percents is attained. At this level, the maximum voltage change potential difference is close to the span of the range of permissible voltages (chosen here to be between 0.95 p.u. and 1.05 p.u.). This observation gives us a hint that we may be able to define a domain of applicability for the criticality criterion through the mere calculation of the maximal difference of active voltage change potential. To investigate this claim, we have completed various sets of simulation while modifying several relevant parameters: the R/X ratio of the lines, the storage discharge time and the number of nodes equipped. For each set of simulation, the level of penetration varies between 10 and 150 percents. The results are presented below. We observe that the number of critical time steps missed remains relatively low (under 60) compared to the total number of time steps (8760) when the maximal difference of active voltage change potential is below the span of permissible voltage range. As the original goal stated is to find the minimal total cost associated with the operation of such a system, not taking into account certain critical time steps will have the consequence to under evaluate this cost. Thus, we compute the difference between the total cost taking into account all the critical time steps and the total cost obtained through the use of the criticality criterion, and define this as the loss of precision incurred by the implementation of the criticality criterion. It is expressed in percentage of the total cost and is presented below as a function of the maximal difference of active voltage change potential, with varying levels of penetration and numbers of nodes equipped. We can observe that the loss of precision that can be attributed to the use of the criticality criterion remains fairly low, even when the maximum difference of active voltage change potential surpasses the span of the permissible voltage range. This can also be used to determine a domain of applicability for the criticality criterion dependent on the level of precision required by the considered application.

t k t k k t k P V P R V , , ,    (19)   t k P K k t k P K k T t P V V V , , max min _ max max        (20)

IV. CONCLUSION

We have presented a methodology aimed at significantly reducing the computational burden associated with the resolution of optimal power flows for distribution systems, in order to render it more practical for planning purposes. We have verified its validity and the computational gains obtained for a practical case study. Afterwards, we have defined an indirect measurement of the domain of validity with a negligible computational cost and have tested it against several extensions of the aforementioned case study. Finally, we have examined the relationship between the loss of precision incurred from the use of the criticality criterion and the position in the domain of validity. Further research has to be done to try and characterize this relationship in order to be able to determine a priori the loss of precision -or an upper bound for it. This will be done by applying this methodology to representative practical distribution networks, with realistic loading conditions and testing several options for the objective function. Another area of improvement lies in the resolution of the slave problem, which could be sped up by drawing on its specificities, especially the fact that it is close to being a quadratically-constrained quadractic programming problem. Finally, this methodology will be integrated in a comprehensive distribution system planning framework adapted to the development of distributed generation and storage

Figure 1 :

 1 Figure 1 : Total active power load for a winter and a summer week

Figure 2 :

 2 Figure 2 : Evolution of the influence of voltage constraints as a function of penetration level

C

  by choosing the minimum value of t C  .

Figure 3 :Figure 4 :

 34 Figure 3 : Criticality criterion ranking, with 17 nodes equipped and a level of penetration of 100%

Figure 5 :

 5 Figure 5 : Computational time spared as a function of the level of penetration

Figure 6 :

 6 Figure 6 : Bus voltages before and after solving the slave problem

Figure 7 :

 7 Figure 7 : Active powers resulting from the slave problem solution

Figure 8 :

 8 Figure 8 : Reactive powers resulting from the slave problem solution

Figure 9 :Figure 10 :

 910 Figure 9 : Gains obtained by the deployment of storage and PV systems, for various levels of penetration and numbers of nodes equipped

Figure 11 :

 11 Figure 11 : Criticality criterion ranking, with 22 nodes equipped and a level of penetration of 140%

Figure 12 :

 12 Figure 12 : Evolution of criticality criterion efficiency and maximum difference of active voltage change potential, as a function of the level of penetration

Figure 13 :

 13 Figure 13 : Number of critical time steps missed, as a function of maximal difference of voltage change potential, for several simulation setups

Figure 14 :

 14 Figure 14 : Loss of precision as a function of the maximal active voltage change potential


Evangéline Pollard
email: epollard@onera.fr

Benjamin Pannetier
email: bpannetier@onera.fr

Michèle Rombaut
email: michele.rombaut@gipsa-lab.inpg.fr

Evangeline Pollard

Convoy detection processing by using the hybrid algorithm (GMCPHD/VS-IMMC-MHT) and dynamic bayesian networks

Keywords: Multitarget Tracking, GMTI, Convoy detection, GMCPHD, VS-IMMC-MHT, dynamic bayesian network

Convoy detection processing by using the hybrid algorithm (GMCPHD/VS-IMMC-MHT) and dynamic bayesian networks

Introduction

In the battlefield surveillance domain, ground target tracking is a first challenging task to assess the situation [START_REF] Hall | An introduction to multisensor data fusion[END_REF]. Data used for tracking comes from Ground Moving Target Indicator (GMTI) sensors which detect moving targets only by measuring their Doppler frequency. The goal is to have a real ground picture: the number of targets, their dynamics, their relationship... in order to discover military events of interests. In this article, we focus on convoy detection. Some studies on convoy detection based on GMTI signatures already exist [START_REF] Corbeil | Data mining of GMTI radar databases[END_REF][START_REF] Klemm | Tracking of convoys by airborne STAP radar[END_REF][START_REF] Koch | Information fusion apects related to GMTI convoy tracking[END_REF], but our purpose is convoy detection by using target tracks. In this context, two steps are proposed : (1) process a hardy multitarget tracking algorithm in order to detect vehicle aggregates with precision in term of cardinality and state estimation, (2) check if the detected aggregates are convoys or not, by introducing other data types (Synthetic Aperture Radar (SAR), video,. . . ) and by using a data fusion method. This purpose is summarized in Figure 1.

Very efficient tracking algorithms exist today and they have to be adaptated to the very complex ground environment. First, the traffic density is very high and generates a large number of measurements. This characteristic eliminates, in our application, Monte-Carlo techniques [START_REF] Kreucher | Multitarget tracking using the joint multitarget probability density[END_REF][START_REF] Angelova | Extended Object Tracking Using Monte Carlo Methods[END_REF][START_REF] Gordon | Bayesian target tracking after group pattern distortion[END_REF]. Moreover measurements are noisy and can contain many false alarms. Also vehicles on the ground are usually quite manoeuvrable over short periods of time according to the sensor scanning time T . Finally, vehicles are detected by the sensor with probability P D and according to the sensor resolution. In other words, when vehicles are very close together, one measurement can be missing generating the spawned targets. This phenomena added to the problem of data association make the classical algorithms, like SDassignment [START_REF] Deb | A generalized S-D assignment algorithm for multisensor-multitarget state estimation[END_REF] or MHT [START_REF] Bar-Shalom | Multitarget-Multisensor Tracking: Principle and Techniques[END_REF] less efficient to track convoys. Mahler by using his work on FInite Set STatistics [START_REF] Mahler | Detecting, tracking and classifying group targets: a unified approach[END_REF] (FISST) and Random Sets. This filter leads to a new class of algorithms [START_REF] Vo | Performance of PHD Based Multi-Target Filters[END_REF] based on the study of joint density probability of the Random Finite Sets (RFS) describing target dynamics and measurements. The first order moment of this RFS, called the intensity function, is the function whose the integral in any region on state space is the expected number of targets in that region. Points with highest density are then expected targets. To improve the number of targets estimation, Mahler proposes a generalization of the PHD called the CPHD [START_REF] Mahler | Multitarget Bayes filtering via firstorder multitarget moments[END_REF], which jointly propagates the intensity function and the entire probability distribution of the number of targets. Under Gaussian assumptions on target dynamics and birth process, Vo proposes a CPHD recursion called the Gaussian Mixture Cardinalized PHD [START_REF] Vo | Analytic Implementations of the Cardinalized Probability Hypothesis Density Filter[END_REF][START_REF] Ulmke | Gaussian Mixture Cardinalized PHD Filter for Ground Moving Target Tracking[END_REF] (GM-CPHD). This approach gives very encouraging results, in particular for the estimation of the number of targets and seems adaptated to convoy tracking.

Nevertheless, as we will show in Table 7, with manoeuvring targets, the GMCPHD has problems with velocity estimation. From this point of view, the GM-CPHD and the IMM-MHT can be seen as complementary algorithms: the first for the estimation of the number of targets and for an approximate position estimation, and the latter can be used to specify state estimation. The proposed hybridization is described in Figure 1. In this approach, we use a special version of the MHT: the VS-IMMC-MHT [START_REF] Pannetier | Multiple Ground Target Tracking with a GMTI Sensor[END_REF] (Variable Structure -Interacting Multiple Model with Constraints -Multiple Hypothesis Tracking) which uses road segment position from Geographical Information System (GIS) to improve the state estimation. Other authors proposed to combine a PHD filter with other filters [START_REF] Lin | Data association combined with the probability hypothesis density filter for multitarget tracking[END_REF]. By using outputs of our algorithms, we are able to detect aggregates with precision. The second step is to define if they are convoys or not.

Before discuss our approach, we define a convoy as a group of vehicles, evolving on the road, having the same dynamics and generally composed of more than two military vehicles. The distance between two vehicles depends on the environment, but most of the time it is over 100m. Giving these restrictions, we want to produce a general convoy model, able to discriminate convoys from a group of vehicles. We use the Dynamic Bayesian Network (DBN) formalism which seems adaptated to this problematic [START_REF] Johansson | Implementation and integration of a Bayesian Network for prediction of tactical intention into a ground target simulator[END_REF].

The paper is organized as follows: Section 2 is a description of the existing GMCPHD filter, Section 3 details how we use this algorithm in a hybrid version, Section 4 explains how DBNs are used for convoy detection. Finally Section 5 describes our simulation and compares results before we conclude in Section 6.

2 Background on the GMCPHD filter

The PHD filter

A Random Finite Set (RFS) is a finite-set valued random variable which can be generally characterized by a discrete probability distribution and a family of joint probability densities representing the existence probabilities of the target set. Considering the RFS of survival targets S k|k-1 between iterations k -1 and k, the RFS of spawned targets B k|k-1 and the RFS of spontaneous birth targets σ k , the global RFS characterizing the multitarget set can be written as:

X k =   ζ∈X k-1 S k|k-1 (ζ)   ∪   ζ∈X k-1 B k|k-1 (ζ)   ∪ σ k (1)
In the same manner, the multitarget set observation Z k can be seen as a global RFS composed by the RFS of measurements originally from the targets X k and by the RFS of false alarms K k :

Z k = x∈X k Θ k (x) ∪ K k (2)
The PHD traditionally evolves in two steps: prediction and estimation that propagate the multitarget posterior density of the target RFS also called the intensity function v. The prediction state is based on the a posteriori intensity function v k-1 at the previous time k -1, the probability P S for a target to survive between times k -1 and k, the transition function f k|k-1 (.|ζ) given the previous state ζ and the intensity of target birth γ k .

v k|k-1 (x) = P s (ζ).f k|k-1 (x|ζ).v k-1 (ζ)dζ + γ k (x)
(3) Knowing the measurement random set Z k , it is possible to update the intensity function as follows:

v k (x) = (1 -P D )v k|k-1 (x)+ z∈Z k P D .g(z|x)v k|k-1 (x) κ k (z) + P D .g(z|ζ)v k|k-1 (ζ)dζ (4)
where g(z|x) is the likelihood of a measurement z knowing the state of a target x, κ k is the clutter intensity which is modeled by a Poisson process.

The GMCPHD filter

The GMCPHD, proposed by Vo [START_REF] Vo | Analytic Implementations of the Cardinalized Probability Hypothesis Density Filter[END_REF], combines a Gaussian mixture model for the intensity function with the Cardinalized generalization of the PHD filter. That means that the posterior target intensity can be written as a Gaussian mixture:

v k (x) = J k i=1 w k,i N (x; m k,i , P k,i ) (5)
where w k,i , m k,i and P k,i are the weight, mean and covariance of the current Gaussians and J k is their number. Moreover, added to the operations ( 3) and ( 4), the probability to have n targets is predicted and estimated in the same way, as, ∀n ∈ N ⋆ ,

p k|k-1 (n) = n j=0 p Γ (n -j)× ∞ l=j C l j P s , v k-1 j 1 -P s , v k-1 l-j 1, v k-1 l p k-1 (l) (6)
with p Γ (nj) the birth probability of (nj) target and C l j the binomial coefficient with parameters (n, j). Following the Bayes theorem, the estimated cardinality distribution p k|k can be written as a likelihood ratio:

p k|k (n) = L(Z k |n) L(Z k ) p k|k-1 (n) (7)
where L(Z k |n) is the likelihood of the measurements set Z k knowing that there are n targets and L(Z k ) is a normalizing constant.

3 The VS-IMMC-MHT / GM-CPHD hybridization

The labeled GMCPHD

In the classical version of the GMCPHD, the problem of track labeling is not considered. Yet, this step is quite important for complex multitarget scenario. Clark and Panta [START_REF] Clark | The GM-PHD Filter Multiple Target Tracker[END_REF][START_REF] Panta | Data Association and Track Management for the Gaussian Mixture Probability Hypothesis Density[END_REF] proposes method but not adaptated to Gaussian mixture and to a large number of targets. Also, we propose, as an alternative, using the track score for the track initialization and in addition to the statistical distance between peak and predicted track to take into account the global weight for the peak to track association.

Let G be the Gaussian set given by the GMCPHD written:

G k = {w k,i , m k,i , P k,i } i∈{1,...,N G k } = G k,1 , . . . , G k,N G k (8) where N G k is the number of Gaussians (N G k > Nk ) at time k.
A track can be defined as a sequence of estimated states describing the dynamics of one target. The goal of tracking is to offer a list of tracks corresponding to all of the targets. That is why this labeling step is necessary in order to provide a track set chosen amongst the Gaussian set G k . A track T k,i is defined at time k by a state xk,i , a covariance P k,i and a score s k,i :

T k,i = {x k,i , P k,i , s k,i } i∈{1,..., Nk } (9)
The track set is finally written:

T k = T k,1 , . . . , T k, Nk (10)
with Nk the estimation of the number of targets given by the GMCPHD. We define a set of association matrices A k of size Nk × N G k to associate the Gaussian set to the tracks. ∀(m, n) ≤ ( Nk , N G k ), an association matrix A k,i is written as:

A k,i (m, n) = 1 if G k,n
can be associated to T k|k-1,m 0 otherwise (11) with T k|k-1,m the predicted track m and knowing that a track is associated at most to one Gaussian. A Gaussian peak n is said associable to a track m if it satisfies a gating test around the predicted position of the track.

We define a weight matrix

W k of size Nk ×N G k defined as follows, ∀(m, n) ≤ ( Nk , N G k ): W k (m, n) = w k,n if G k,n can be ass. to T k|k-1,m 0 otherwise (12) If Nk > Nk-1 ,
one or more new tracks must be initialized and each Gaussian is a potential new track. In matrix

W k , ∀m ∈ {1, . . . , N G k }, ∀l ∈ { Nk-1 +1, . . . , Nk }, W k (m, l) = w k,l (13)
In the same way, if Nk < Nk-1 , some tracks must be deleted. Weakly weighted tracks cannot be deleted because of the detection probability, which is why tracks with the lowest score are deleted.

Finally, we compute the set of global weight of an association:

W g k = Nk m=1 N G k n=1 A k (m, n).W k (m, n) (14)
And the association matrices which maximize the weight are written as:

A ⋆ k = argmax A k W g k (15)
Similarly, the cost matrix

C k of size Nk × N G k is writ- ten as, ∀(m, n) \ ( Nk , N G k ), C k (m, n) = c(m, n) if G k,n
can be ass. to T k|k-1,m 0 otherwise (16) with c(m, n) the cost of the association of the predicted track m with the Gaussian n written as the negative Napierian logarithm of the likelihood ratio,

∀(m, n) ∈ ( Nk , N G k ), c(m, n) = -ln P D .Λ(G k,n |T k|k-1,m ) β F A (17)
with β F A the spatial false alarm density and Λ(G k,n ) the likelihood of the Gaussian n knowing the predicted position of the track m, calculated as a Gaussian density.

Finally, the global association cost is computed as:

C g k = Nk m=1 N G k n=1 A ⋆ k (m, n).C k (m, n) (18)
And the best association A ⋆⋆ is computed like the minimal cost matrix:

A ⋆⋆ k = argmin A ⋆ k C g k (19)

The hybridization

The GMCPHD produces a reliably estimation of the number of targets, whereas the VS-IMMC-MHT is effective to give a good estimation of the target state by introducing road coordinates when targets are not close together, because of the problem for MHT algorithm to evaluate the number of targets. We propose therefore to use these two algorithms as complementary filters: the first estimates the number of targets and the approximate target position and the second increases the accuracy for the target state estimation. The two algorithms are running simultaneously. Then, a gating process is applied around the target position given by the GMCPHD, to select MHT tracks. Finally, MHT tracks which have the highest score are selected. If a PHD track is not associated to any MHT track, the GMCPHD track is kept.

This approach combines the advantages of the different algorithms without increasing the processing time:

• Robust to target maneuvers by using IMM

• Good precision for state estimation by using road coordinates

• Good estimation of the number of targets

• No performance decrease when targets are close together Different algorithms performances are compared in Section 5 in a complex scenario. But before let us define the proposed convoy detection method.

Description of a convoy 4.1 Some definitions

A convoy is defined as a vehicle set evolving approximatively with the same dynamics during a long time. These vehicles are moving on the road under a limited velocity (<20m/s). They must stay at sight with almost constant distances between them (mostly 100m). Criteria describing a convoy are manifold and of different natures, moreover variables are discrete. That is why, bayesian networks represent an interesting formalism in our application as in similar thematics [START_REF] Okello | Threat assessment using Bayesian networks[END_REF][START_REF] Denis | Spatio-temporal pattern detection using dynamic Bayesian networks[END_REF][START_REF] Johansson | Implementation and integration of a Bayesian Network for prediction of tactical intention into a ground target simulator[END_REF][START_REF] Singh | Modeling threats[END_REF].

A Bayesian Network (BN) is a graphical model for representing dependency relation between a set of random variables. Graphically, each variable is represented by a node and an arc, from a node X i to a node X j , means that X i "causes" X j , ∀(i, j) ∈ {1, . . . , N } 2 . Finally, the joint probability is computed as:

P (X 1 , . . . , X N ) = N i=1 P (X i |P a(X i )) (20)
where P a(X i ) are parent nodes of X i . The Dynamic Bayesian Networks (DBN) are an extension of BN, which take into account the time evolution of random variables. The convoy detection approach is bounded to the time evolution as shown in Figure 2. For example, variable X 5 is time depending, because the type information can come from heterogeneous sources (SAR, video, . . . ) with different scanning times, and variable X 9 is confirmed with time.

X 1 X 2 X 3 X 4 X k 5 X 8 X 6 X 7 X k 9 X1: Velocity < 80km/h

Conditional Probability Distribution evaluation

If independency relation between variables can be very intuitively established, one difficulty with DBN is to evaluate the Conditional Probability Distribution (CPD) of each node given its parents. If data sets are available, these prior probabilities can be learned, but in our case, they are evaluated by experts, according to a certain weight to each parameters. For example, if a convoy is detected at time k -1, the probability to detect one at time k is high and the prior probability given to this parameter must be "relatively" high. As described in [START_REF] Benavoli | An approach to threat assessment based on evidential networks[END_REF], we propose heuristic rules to represent relationships between variables :

5 × X k-1 9 + 1.5 × X 3 + X 4 + X 5 + 1.5 × X 8 = X k 9 (21)
As said, this rule means that the probability to have a convoy at time k is half-depending on the fact to have a convoy at time k -1, and that we care more criteria on distance and velocity than criteria "on road" or "vehicle type".

Probability transformation

Another difficulty is the transformation of numerical data (number of targets, target position and velocity, road position) into a probability. This step is done by using probability distributions or fuzzy transformations like linear transformation (cf. • p(X 1 ) is computed according to a Rayleigh distribution.

• p(X 2 ) is following a fuzzy linear transform using the difference between velocity mean at time k and at previous times.

• p(X 4 ) is computed according to a χ 2 distribution.

• p(X 6 ) is computed using a fuzzy gaussian transformation by studding the distribution of distances between vehicles of the aggregate.

• p(X 7 ) is computed using a fuzzy linear transformation by examining the variation of the convoy length over time

Inference

The next step consists to propagate the information through the network. It is called the inference. Many algorithms exist like JLO [START_REF] Jensen | Bayesian updating in recursive graphical models by local computation[END_REF] from the names of its authors or Expectation-Maximization (EM) algorithm. We choose arbitrarily the JLO algorithm adaptated to discrete nodes and available in the Murphy's Bayes net toolbox [START_REF] Murphy | Bayes net toolbox[END_REF].

Targets number estimation

Computing the probability p(X 9 ) for an aggregate to be a convoy is a first step (cf. Figure 9), but it is possible to take into account the average number of targets belonging to the convoy. First, we know the number of targets in the aggregate and moving in the same direction. If for instance at time k, we detect N (k) = 5, while there was 4 until there, we have to propagate the information and to compute simultaneously the probability to have a convoy with 5 vehicles, and a convoy with the 4 best located target tracks. Mathematically, it means we compute p(X 9 , N C ), with N C the set of different values taken by N k , where N k = {N (1), . . . , N (k)} is the sequence of mean number of targets in the aggregate, moving in the same direction.

However, as shown in Figure 11, it is not easy to discriminate certain cases, here the cases N C = 5 and N C = 6 (the reality is N C = 6). If, at the beginning of the simulation, we detect N C = 5, we must continue to compute the probability to have a 5 target convoy, because we are possibly in the case of an overtaking, but it is not realistic, to support this assumption against the 6 target convoy if the sequence of measurement never gives again N (k) = 5. That is why we introduce the local estimated cardinality of the Gaussian mixture on the aggregate surface, computed as

N C k = Nmax i=1
i.p C k|k (i) knowing the sequence of average number of targets.

Finally, the probability becomes:

p(X k 9 , N C , N C k |N k ) = p(X 9 , N C ).p( N C k |N k ) (22)
By considering a Markovian assumption and Bayes theorem, the probability is computed as:

p( N C k |N k ) = p(N (k)| N C k , N (k -1)).p( N C k |N (k -1)). 1 c (23) with c a normalization constant, p( N C k |N k-1 ) = N ( N C k ; N k-1 , σ 2 N
) is computed as the normal density with mean N k-1 and variance σ 2 N , and

p(N (k)| N C k , N k-1
) is computed by using a linear transformation.

Simulation and results

In the following, we present some simulation results that illustrate the performances of the proposed hybridization. These are compared to the performances of a classical IMM-MHT, a labeled GMCPHD, a VS-IMMC-MHT and an hybridization GMCPHD/IMM-MHT. Then we present some results on the convoy detection.

Scenario

The GMTI sensor has a linear trajectory, its velocity is 30m/s and its altitude is 4000m. The typical measurement error is 20m in range and 0.008rad in azimuth. The sensor scan time is T = 10s. Scenario time is limited to 500s. The false alarm density is β F A = 8.92.10 -9 and the detection probability is P D = 0.9. Target trajectories are illustrated in Figure 4, while cumulated MTI reports are shown in Figure 5. In the scenario, one 6 target convoy (Target 1-6) is moving on the main road with a constant velocity of 10m/s from South to North. An independant target (Target 7) is moving on the same road in the same direction but with a constant velocity of 15m/s and overtakes the convoy between time t=150s and t=350s approximately.

Results

The performances of tracking algorithms have been compared for 100 independent Monte Carlo runs. Fig- Concerning the convoy probability, p(X 9 ) is evolving progressively from 0.5 to 0.6 (cf. Figure 9) with some picks which indicate a change of cardinality in the aggregate. By introducing N C , we begin to estimate the number of targets in the aggregate (cf. Figure 11). We discriminate the case N C = 7, but we cannot decide between N C = 5 or 6. Finally, by introducing N C k knowing N k (cf. Figure 10), the case N C = 6 appears as the

Conclusion

The new approach for convoy detection has shown its efficiency on a complex multitarget scenario. Several theoretical contributions have been proposed. The first one concerns the labeled version of the GMCPHD that allows to differentiate the tracks. The second con- tribution concerns the hybridization of the GMCPHD algorithm to the VS-IMMC-MHT algorithm in order to improve the performances, specially for group of closely spaced objects. Finally, the third contribution concerns the convoy model by using DBN that proposes an original answer to convoy detection process. This has been tested on several scenarios not presented in the paper. The next step is now the problem of make a decision which stays entire.

Figure 1 :

 1 Figure 1: Convoy detection process

Figure 2 :

 2 Figure 2: Dynamic bayesian network for convoy detection. The gray nodes represent states depending on their previous state.

Figure 3 :

 3 Figure 3: Examples of transformation

  Fig 3 (a)) or Gaussian transformation (cf. Fig 3 (b)):

Figure 4 :

 4 Figure 4: Scenario

Figure 5 :

 5 Figure 5: Cumulated MTI reports The simulation parameters are presented in Tables 1 to 5. Name Value CV model noise 1 0.05m.s -2 CV model noise 2 0.8m.s -2 Model noise 3 (STOP) 0m.s -2

ure 6

 6 shows the average RMSE (Root Mean Square Error) of each target in position, Figure7average RMSE in velocity and Figure8is the track length ratio of each target. The IMM-MHT offers acceptable performances in state estimation, while the VS-IMMC-MHT improves highly position estimation. The GMCPHD produces lower performances in term of state estimation, but the track length ratio is close to 1. The hybrid version (Hybrid 1 is the hybridization of IMM-MHT and GMCPHD, Hybrid 2 is the hybridization of VS-IMMC-IMM-MHT and GMCPHD) is a good compromise between the two sorts of algorithms. The track length ratios have similar values as the GMCPHD, whereas, the state estimation are similar for Hybrid 1 to the IMM-MHT and for Hybrid 2 to the VS-IMMC-MHT.

Figure 6 :

 6 Figure 6: Average RMSE in position of each target

Figure 7 :

 7 Figure 7: Average RMSE in velocity of each target

Figure 8 :

 8 Figure 8: Track length ratio

Figure 9 :

 9 Figure 9: p(X 9 )

Table 1 :

 1 The IMM parameters

Table 2 :

 2 The MHT parameters

	Name	Value
	Survival probability	0.98
	Initial Gaussian weight	10 -3
	Pruning threshold	10 -2
	Merging threshold	20
	Maximum number of targets	50
	Maximum number of Gaussians	50
	Average number of birth
	Model noise	2
	Maximum velocity	20

Table 3 :

 3 The GMCPHD parameters

	Name	Value
	CV model noise 1 in normal direction	0.1
	CV model noise 1 in orthogonal direction	0.1
	CV model noise 2 in normal direction	0.6
	CV model noise 2 in orthogonal direction	0.4
	Maximum value for off road velocity	9m.s -1

Table 4 :

 4 The VS-IMMC parameters

	Name	Value
	Number of iterations for score calculation	3
	Weight threshold for new track	0.8

Table 5 :

 5 The hybridization parameters