Abdallah Zegaoui
email: abd.zegaoui@gmail.com

Mohamed Boutoubat

Jean-Paul Sawicki

Fatma Zohra Kessaissia

Abdelkader Djahbar

Michel Aillerie
email: aillerie@metz.supelec.fr

Enhanced Model of Photovoltaic Cell/Panel/Array Considering the Direct and Reverse Modes

Keywords: Photovoltaic cells, Photovoltaic panels, PV cell/panel modeling, Direct and reverse operating mode, I-V characteristics, Environmental parameters

This paper presents an improved generalized physical model for photovoltaic, PV cells, panels and arrays taking into account the behavior of these devices when considering their biasing existing in direct and reverse modes. Existing PV physical models generally are very efficient for simulating influence of irradiation changes on the short circuit current but they could not visualize the influences of temperature changes. The Enhanced Direct and Reverse Mode model, named EDRM model, enlightens the influence on the short-circuit current of both temperature and irradiation in the reverse mode of the considered PV devices. Due to its easy implementation, the proposed model can be a useful power tool for the development of new photovoltaic systems taking into account and in a more exhaustive manner, environmental conditions. The developed model was tested on a marketed PV panel and it gives a satisfactory results compared with parameters given in the manufacturer datasheet.

INTRODUCTION AND INITIAL BACKGROUND

The efficiency of photovoltaic (PV) panels and arrays primary depends on the efficiency of the basic elements being the PV cells. A PV cell can operate in the first quadrant of the I-V characteristics as generator or in the second one as receiver. When on PV cell in a string work in this last operating mode, it imposes a reverse voltage against the others PV cells of the string. As shown in Ref. [START_REF] Naka | Simulation of i-v characteristics of a pv module with shaded pv cells[END_REF][START_REF] Tu | Development of generalized photovoltaic model using matlab/Simulink[END_REF][START_REF] Bishop | Computer simulation of the effects of electrical mismatches in photovoltaic cell interconnection circuit[END_REF][START_REF] Campbell | A circuit-based photovoltaic array model for power studies[END_REF][START_REF] Gow | Development of photovoltaic array model for use in power-electronics simulation studies[END_REF][START_REF] Zegaoui | Simulation of photovoltaic generators and comparison of two common maximum power point trackers[END_REF] and by Zaraket et al. in Ref. [START_REF] Zaraket | The Effect of Electrical stress under temperature in the characteristics of PV Solar Modules[END_REF][START_REF] Zaraket | Dark and illuminated characteristics of photovoltaic solar modules. Part I: Influence of dark electrical stress[END_REF][START_REF] Zaraket | Dark and illuminated characteristics of photovoltaic solar modules. Part II: Influence of light electrical stress[END_REF][START_REF] Zaraket | Parameters and characteristics of PV solar modules under the influence of thermal stresses[END_REF] who have experimentally studying the influence of electrical and thermal stresses on cell efficiency, shadow induces the PV cell to act as load and to consume the power delivered from other un-shading cells, thus inducing a rising of its temperature reaching creation of hot spot. Therefore, the authors in [START_REF] Zegaoui | Photovoltaic Cell/Panel/Array Characterizations and Modelling Considering both Reverse and Direct Modes[END_REF] have originally proposed a direct and reverse mode model (named DRM model) able to describe I-V characteristics in both direct and reverse operating mode according to the experiments made on the reverse mode biasing of 60 polycrystalline silicon cells.

The existing models, except the DRM model but even the relatively recent ones based on transistor instead of diode [START_REF] Zegaoui | Universal Transistor-based hardware SIMulator for real time simulation of photovoltaic generators[END_REF][START_REF] Petit | The Transistor based Direct and Reverse Mode model for photovoltaic strings and panels[END_REF], do not take into account the reverse operating mode. DRM model accurately reproduce the I-V behavior of PV cells, PV panel and PV array according to the irradiation and the temperature changes. Nevertheless, the influence of the temperature variation is not visible on the I-V characteristics and the influence of temperature changes on the short current in the second operating quadrant was not considered avoiding the consideration of the environmental parameters of the cell in operation.

In the current work and contribution, we present an enhanced DRM model, (EDRM model) for PV cells/panel/array modelling named EDRM considering the influence of temperature variations. The enhanced model hugely improve the existing one as it can be apply to cells or panel or array, not characterized in static mode in laboratory but characterized in real operation modes.

The development done for the EDRM model was supported by experimental I-V characterizations of 60 polycrystalline PV cells in reverse mode produced from the same initial ingot. The I-V characteristics were shown in Fig . 1. Characteristics extracted from Ref. [START_REF] Zegaoui | Photovoltaic Cell/Panel/Array Characterizations and Modelling Considering both Reverse and Direct Modes[END_REF] From figure 1, we note the presence of a huge discrepancy in the I-V behavior between all PV cells even if the 60 cells provide from the same ingot. We note also the appearance of two avalanche voltages on the reverse mode: the first value is around -5V and the second value close to -14V. These experimental finding had brought to propose in Ref. [START_REF] Zegaoui | Photovoltaic Cell/Panel/Array Characterizations and Modelling Considering both Reverse and Direct Modes[END_REF] the initial basic DRM model as reported hereafter in Fig. 2 to obtain a self-consistent publication. This model was developed under Orcad/Pspice as shown in Fig. 3.

FIGURE 2.

The basic DRM model originally developed in Ref. [START_REF] Zegaoui | Photovoltaic Cell/Panel/Array Characterizations and Modelling Considering both Reverse and Direct Modes[END_REF].

FIGURE 3.

The basic DRM model implemented under Orcad/Pspice (from Ref. [START_REF] Zegaoui | Photovoltaic Cell/Panel/Array Characterizations and Modelling Considering both Reverse and Direct Modes[END_REF]). To test the effectiveness of our approach when temperature changes are considered in the behavior of PV elements, the EDRM model was applied to a PV cell, with the following steps of temperature: -10, 0, 25, 50, 70, 85 and 100 °C. Fig. 5 shows the I-V characteristics of this PV cell as function of temperature. In Fig. 5, one can remark that the rise in temperature causes a decrease of the open circuit voltage inducing a decrease of the maximum power of the PV, as expected by the usual behavior of PV cell with temperature. The the rise in temperature is also at the origin of an increase of the short-circuit current, which directly correspond to an increase of the photocurrent with the temperature. For the specific simulated PV cell, in the temperature range, the short-circuit current rises from 1.5 to up than 1.6 amps for the entire region that the PV cell behaves as a generator, i.e. when it remains in the first quadrant on the I-V characteristics. In the reverse mode biasing, the current keeps the short current values depending to the temperature changes. Nevertheless, the reverse voltage imposed by the PV cell in this second quadrant keeps going up.

THE NEW ENHANCED DRM MODEL

THE EDRM MODEL OF A SHADED MARKETED PV PANEL

In order to show the influence of both shadow and temperature changes, we have shaded the entire effective surface of one PV cell in the PV panel and then varied the temperature steps by steps as simulated in Fig. 5. The corresponding EDRM model of the simulated PV panel was implemented under Orcad/Pspice software. It corresponds to an extension of the EDRM model shown in Fig. 4 for one cell. Indeed, for the modeling of a panel, the model is simply built with an association in series and in parallel of the EDRM model for one cell according to the datasheet of the solar module configuration under test.

The BP solar module BP350J is a 50Wc PV panel composed by two parallel strings of 19 polycrystalline silicon PV cells. Its open voltage value and short current value are 21.8 V and 3 A under STD operating conditions, respectively. The maximum power point coordinates in the I-V characteristics plane are (2.9 A, 17.5 V). Only one by-pass diode is mounted for its protection against inverse current that could be induced by hot spots due to shading on the surface of the panel. The I-V characteristics, under temperature variations, as provided by manufacturer in the datasheet, are shown in Fig. 6. Simulated results of I-V and P-V characteristic curves due to the temperature variations were presented in Figs. 7 and8, respectively. Mainly in Fig. 7, we present the influence of the temperature variations on the I-V characteristics of the simulated shadowed BP350J PV BP solar with a zoom around the short circuit current area; whereas in Fig. 8, we present the displacement of the local MPP on the P-V characteristics of the BP350J solar panel under partial shadow and temperature changes. Some pertinent information can be extracted of these figures. At first, as it is well-known, and independently of the temperature values, these figures enlighten the direct effect of the shadow on the global energetic performance of the panel. Indeed, the electrical characteristics present knees for the I-V ones and power local maximums on the P-V ones. It is to be noted that the presence of local maximum power points is at the origin of large problem in the operation of converters or inverters whose one main function is to track the maximum power point so as to convert the maximum power. At second, these figures confirm, as it is also well-known, that the temperature rise induce a global decrease of the delivered power. Finally, the efficiency of the EDRM model and it pertinence in the modeling of shading effect are confirmed by the change with temperature of the short circuit current and on the open voltage circuit, observed in Figs. 7 and8 with a faithful reproduction of values given by the manufacturer and shown in Fig. 6. operate in the reverse biased mode and experimental results are compared to the simulated ones. The EDRM model is able to reproduce the influence of environmental variables, i.e. temperature and irradiation on the I-V and P-V characteristics in both the first and the second quadrant of the operating modes. The are affected by temperature variations, our simulations are validated by the fact that the changes of the photocurrent and as a result the short circuit current are in the same direction as the temperature variations. To prove the validity of the proposed model, the degree of the temperature influence on the photocurrent is compared to that mentioned on the I-V curves of a marketed PV panel. We have shown that the EDRM model faithfully reproduce the characteristics expected in the case of shading with the values of the short-circuit current and the open circuit voltage given by the manufacturer. An extension in applicative possibilities of the EDRM model is the modeling of PV arrays, which is easily possible using a series-parallel combination of the EDRM model of a panel.

FIGURE 1 .

 1 FIGURE 1. Experimental I-V curves of 60 polycrystalline PV cells in the reverse mode biasing.Characteristics extracted from Ref.[START_REF] Zegaoui | Photovoltaic Cell/Panel/Array Characterizations and Modelling Considering both Reverse and Direct Modes[END_REF]

Figure 4

 4 Figure 4 show the electronic diagram of the new enhanced DRM model of PV cells. This model improve the basic DRM model as it takes into consideration in the reverse operating mode both the variation of the

  photocurrent and the short circuit current due to environmental parameters' changes, especially link to the influence of temperature variations. The main modification of the basic DRM model is the adding of a diode (D384) able to perfectly reproduce the influence of temperature variations on the photocurrent on the I-V characteristics and on the short current.

FIGURE 4 .

 4 FIGURE 4. The EDRM model for PV cell implemented under Orcad/Pspice.

FIGURE 5 .

 5 FIGURE 5. I-V characteristics of a PV cell in both direct and reverse mode biasing under influence of temperature changes by applying the EDRM model.

FIGURE 6 .

 6 FIGURE 6. Influence of temperature variations on the I-V characteristics of PV BP solar BP350J from manufacturer datasheet.

Technologies and Materials for Renewable Energy, Environment and Sustainability AIP Conf. Proc. 1968, 030015-1-030015-6; https://doi.org/10.1063/1.5039202 Published by AIP Publishing. 978-0-7354-1675-8/$30.00

The temperature influence on the current is evaluated thanks to the EDRM model by a value equal to 0.06445 mA/°C. So, the short circuit current is affected by a rate value equal to 0.064%/°C. This value is nearby the value mentioned by the manufacturer in the panel datasheet, which is 0.065%/°C. This result proves the faithfulness of the new proposed EDRM model for modeling PV cells, but also, as shown here, PV panels obtained by duplicating the model for one cell and one can easily imagine that the EDRM model will be also faithful for the modeling of PV arrays extending the EDRM model presented above for a panel, with now, the respect of the specific seriesparallel configuration of the wide array.

CONCLUSION

In this work, a new general enhanced physical model for PV cells, panels and arrays operating in both direct and reverse biasing modes is proposed. The new model, named EDRM, is based on the use of a conventional one or two diode models used for modeling of PV cell. Experiments were performed on polycrystalline PV cells


Guillaume Carlier
email: carlier@ceremade.dauphine.fr

Teresa Radice
email: teresa.radice@unina.it

Approximation of variational problems with a convexity constraint by PDEs of Abreu type

Keywords: Abreu equation, Monge-Ampère operator, calculus of variations with a convexity constraint. MS Classification: 35G30, 49K30

Motivated by some variational problems subject to a convexity constraint, we consider an approximation using the logarithm of the Hessian determinant as a barrier for the constraint. We show that the minimizer of this penalization can be approached by solving a second boundary value problem for Abreu's equation which is a well-posed nonlinear fourth-order elliptic problem. More interestingly, a similar approximation result holds for the initial constrained variational problem.

Introduction

Given Ω, a bounded, open, convex subset of R d with d ≥ 2, F : Ω × R → R strictly convex in its second argument, and ϕ a uniformly convex and smooth function defined in a neighbourhood of Ω, we are interested in the variational problem with a convexity constraint:

inf u∈S[ϕ,Ω] J 0 (u) := Ω F (x, u(x))dx (1.1)
where S[ϕ, Ω] consists of all convex functions on Ω which admit a convex extension by ϕ in a neighbourhood of Ω. This is a way to express in some weak sense the boundary conditions

u = ϕ and ∂ ν u ≤ ∂ ν ϕ on ∂Ω, (1.2)
where ν denotes the outward normal to ∂Ω and ∂ ν denotes the normal derivative.

Due to the convexity constraint, it is really difficult to write a tractable Euler-Lagrange equation for (1.1) (see [START_REF] Lions | Identification du cône dual des fonctions convexes et applications[END_REF], [START_REF] Carlier | Representation of the polar cone of convex functions and applications[END_REF]). One may therefore wish to construct suitable penalizations for the convexity constraint which force the minimizers to somehow remain in the interior of the constraint and thus to be a critical point of the penalized functional. Since the seminal work of Trudinger and Wang [START_REF] Trudinger | The affine Plateau problem[END_REF][START_REF] Trudinger | Boundary regularity for the Monge-Ampère and affine maximal surface equations[END_REF] on the prescribed affine mean curvature equation, the regularity of convex solutions of fourth-order nonlinear PDEs which are Euler-Lagrange equations of convex functionals involving the Hessian determinant have received a lot of attention. In particular, the Abreu equation which corresponds to the logarithm of the Hessian determinant has been studied by Zhou [START_REF] Zhou | The first boundary value problem for Abreu's equation[END_REF] in dimension 2 and more recently by Chau and Weinkove [START_REF] Chau | Monge-Ampère functionals and the second boundary value problem[END_REF] and Le [START_REF] Le | Global second derivative estimates for the second boundary value problem of the prescribed affine mean curvature and Abreu's equations[END_REF][START_REF] Le | W 4,p solution to the second boundary value problem of the prescribed affine mean curvature and Abreu's equations[END_REF] in higher dimensions. What the well-posedness and regularity results of these references in particular suggest is that a penalization involving the logarithm of the Hessian determinant should act as a good barrier for the convexity constraint in problems like (1.1). This was indeed confirmed numerically at a discretized level, see [START_REF] Benamou | Discretization of functionals involving the Monge-Ampère operator[END_REF].

Our goal is precisely to show that one can indeed approximate (1.1) by a suitable boundary value problem for the Abreu equation. To do so, we first introduce a penalized version of (1.1) with a small parameter ε > 0:

inf v∈S[ϕ,Ω] J ε (v) := J 0 (v) -εF Ω (v) (1.3)
where, when v ∈ S[ϕ, Ω] is smooth and strongly convex, (see section 2 for the definition for an arbitrary v ∈ S[ϕ, Ω]), F Ω (v) is defined by

F Ω (v) := Ω log(det D 2 v).
Using the convexity of J ε setting f (x, u) := ∂ u F (x, u), one can easily see that if u is smooth and uniformly convex up to ∂Ω, and solves the first-boundary problem for Abreu equation

εU ij w ij = f (x, u) in Ω, u = ϕ and ∂ ν u = ∂ ν ϕ on ∂Ω (1.4)
where w := det(D 2 u) -1 and U denotes the cofactor matrix of D 2 u then it is indeed the solution of (1.3). It turns out however that the second-boundary value problem (where instead of prescribing both values of u and ∇u one rather prescribes u and det(D 2 u) on ∂Ω) is much more well-behaved, see [START_REF] Chau | Monge-Ampère functionals and the second boundary value problem[END_REF][START_REF] Le | Global second derivative estimates for the second boundary value problem of the prescribed affine mean curvature and Abreu's equations[END_REF][START_REF] Le | W 4,p solution to the second boundary value problem of the prescribed affine mean curvature and Abreu's equations[END_REF] and it was indeed used as an approximation for the affine Plateau problem in [START_REF] Trudinger | The affine Plateau problem[END_REF]. We shall also consider an extra approximation parameter and a second-boundary value problem on a larger domain and show that it approximates correctly not only (1.3) but also the intial problem (1.1) as the parameter converges to zero. The paper is organized as follows. Section 2 gives some preliminaries. In section 3, we show a Γ-convergence result for J ε . In section 4, we consider an approximation by a second boundary value problem on a ball B containing Ω, with a further penalization 1 δ (u-ϕ) on B\Ω, for which we prove existence and uniqueness of a smooth solution. In section 5, we show that when δ → 0, we recover the minimizer of the problem from section 3. Finally, we also show full convergence of the second boundary value problem to the initial constrained variational problem (1.1) when δ = δ ε → 0 as ε → 0, provided F satisfies a suitable uniform convexity assumption.

Preliminaries

In the sequel, Ω will be an open, bounded and convex subset of R d , d ≥ 2. We are also given an open ball B containing Ω and assume that the boundary datum ϕ satisfies for some λ > 0:

ϕ ∈ C 3,1 (B), ϕ = 0 on ∂B, D 2 ϕ ≥ λ id on B.
(2.1)

We then define S[ϕ, Ω] as the set of convex functions on Ω, which, once extended by ϕ on B \ Ω, are convex on B. Note that elements of S[ϕ, Ω] coincide with ϕ on ∂Ω and are Lipschitz continuous with Lipschitz constant at most ∇ϕ L ∞ (B) so that S[ϕ, Ω] is compact for the topology of uniform convergence.

Finally, we assume that the integrand F : (x, u) ∈ Ω × R → F (x, u) in the definition of J 0 in (1.1) is measurable with respect to x, strictly convex and differentiable with respect to u and such that that F (., 0) ∈ L 1 (Ω) and f (x, u) := ∂ u F (x, u) satisfies f (., u) ∈ L ∞ (Ω) for every u ∈ R. These assumptions in particular guarantee that the convex functional J 0 is everywhere continuous and Gâteaux differentiable on S[ϕ, Ω].

Following [START_REF] Trudinger | The affine Plateau problem[END_REF][START_REF] Trudinger | Boundary regularity for the Monge-Ampère and affine maximal surface equations[END_REF][START_REF] Zhou | The first boundary value problem for Abreu's equation[END_REF], let us recall how to define F Ω (v) for an arbitrary convex function v on Ω, first recall that the subdifferential of v at x ∈ Ω is given by

∂v(x) := {p ∈ R d : v(y) -v(x) ≥ p • (y -x), ∀y ∈ Ω}.
The Monge-Ampère measure of v, denoted µ[v] is then defined by

µ[v](A) := |∂v(A)|
for every Borel subset A of Ω. From the seminal results of Alexandrov (see [START_REF] Cristian | The Monge-Ampère equation[END_REF]), µ[v] is indeed a Radon measure and v → µ[v] is weakly continuous in the sense that whenever v n are convex functions which locally uniformly converge to v then lim sup

n µ[v n ](F ) ≤ µ[v](F ), ∀ F ⊂ Ω, closed .
Decomposing the Monge-Ampère measure into its absolutely continuous part and its singular part (with respect to the Lebesgue measure L d ) as

µ[v] = µ r [v] + µ s [v], µ r [v] ≪ L d , µ s [v] ⊥ L d .
Thanks to Alexandrov's theorem, v is differentiable twice a.e., at such points of twice differentiability, we denote by ∂ 2 v its Hessian matrix, Trudinger and Wang proved in [START_REF] Trudinger | The affine Plateau problem[END_REF] that det(∂ 2 v) is the density of µ r [v] with respect to L d , and following their approach, one can define the functional F Ω by

F Ω (v) := Ω log(det ∂ 2 v(x))dx, ∀v ∈ S[ϕ, Ω]. (2.2)
It is well-known that F Ω is a concave functional and we refer to [START_REF] Trudinger | The Bernstein problem for affine maximal hypersurfaces[END_REF][START_REF] Trudinger | The affine Plateau problem[END_REF][START_REF] Zhou | The first boundary value problem for Abreu's equation[END_REF] for a proof of the useful properties of F Ω recalled below in Lemmas 2.1 and 2.2

Lemma 2.1. The functional v ∈ S[ϕ, Ω] → F Ω (v) defined in (2.
2) is concave, upper semi-continuous for the topology of local uniform convergence and bounded from above on S[ϕ, Ω] with the explicit bound (where c d denotes the measure of the unit ball of R d )

F Ω (v) ≤ C Ω,ϕ := |Ω| log c d ∇ϕ d L ∞ |Ω| , ∀v ∈ S[ϕ, Ω]. (2.3)
As we shall also work on the larger domain B, it will be also convenient to consider for every open subset ω of B and every convex function u on B the concave functional

F ω (v) := ω log(det ∂ 2 v(x))dx.
(2.4)

Following the same lines as Lemma 6.4 in Trudinger-Wang [START_REF] Trudinger | The Bernstein problem for affine maximal hypersurfaces[END_REF], we also have:

F ω (u n ) ≤ F ω (u).
3 Logarithmic penalization

Given ε > 0, we consider inf v∈S[ϕ,Ω] J ε (v) := J 0 (v) -εF Ω (v). (3.1)
Since J ε is strictly convex and lsc on the convex compact set S[ϕ, Ω], we immediately have:

Proposition 3.1. Problem (3.1) admits a unique solution v ε .
Arguing exactly as in [START_REF] Trudinger | The affine Plateau problem[END_REF][START_REF] Zhou | The first boundary value problem for Abreu's equation[END_REF] by using Alexandrov's maximum principle, one can show: Lemma 3.2. Let ε > 0 and v ε be the solution of (3.1) then µ s [v ε ] = 0 i.e. µ[v ε ] has no singular part. Remark 3.3. Let us remark that Lemma 3.2 enables one to express -F Ω (v ε ) in an alternative way as the entropy of the push-forward of the Lebesgue measure on Ω by ∇v ε . Also, thanks to Lemma 3.2, one can prove uniqueness of the solution of (3.1) when J 0 is convex but not necessarily strictly convex.

In dimension 2, we actually even have a uniform local bound on det(∂ 2 v ε ):

Proposition 3.4. Let d = 2, ε > 0 and v ε be the solution of (3.1), then µ[v ε ] = det(∂ 2 v ε ) ∈ L ∞ loc (Ω).
Proof. It follows from Theorem 5.1 and Proposition 4.3 that v ε is the uniform limit as δ → 0 of a sequence of smooth functions

(v δ ε ) in S[ϕ, Ω] such that, for every open subset ω with ω ⊂⊂ Ω, det(D 2 v δ ε ) L ∞ (ω) ≤ C
where C is a constant that depends on ε and ω but not on δ . By weak convergence of Monge-Ampère measures we deduce that det(∂

2 v ε ) ∈ L ∞ loc (Ω).
Let us now state a Γ-convergence result for J ε :

Proposition 3.5. The family of functionals J ε defined on S[ϕ, Ω] equipped with the topology of uniform convergence Γ-converges to J 0 in particular v ε converges uniformly to the solution of (1.1).

Proof. Assume u ε is a family in S[ϕ, Ω] that converges uniformly as ε → 0 to u, thanks to (2.3) and Fatou's Lemma, we have

lim inf ε J ε (u ε ) ≥ lim inf ε (J 0 (u ε ) -εC Ω,ϕ ) ≥ J 0 (u).
Given u ∈ S[ϕ, Ω], we now look for a recovery sequence u ε ∈ S[ϕ, Ω] converging to u and such that lim sup ε J ε (u ε ) ≤ J 0 (u), we simply take

u ε := (1 -ε)u + εϕ since ∂ 2 u ε ≥ εD 2 ϕ we have F Ω (u ε ) ≥ d|Ω| log(ε) + Ω log(det(D 2 ϕ))
with the convexity of J 0 , we then have lim sup

ε J ε (u ε ) ≤ lim sup ε ((1 -ε)J 0 (u) + εJ 0 (ϕ)) + O(ε log(ε))) ≤ J 0 (u).

Second boundary value approximation

Having Proposition 3.5 in mind, we now fix the value of ε. Throughout this section, to simplify notations, we therefore take ε = 1 and we are interested in approximating the solution of inf

v∈S[ϕ,Ω] J 1 (v) := Ω F (x, v(x))dx -F Ω (v), (4.1)
by a second-boundary value problem for Abreu equation. More precisely given δ > 0, we consider

U ij w ij = f δ (x, u), in B, u = ϕ, w = ψ on ∂B (4.2)
where ψ := det((D 2 ϕ) -1 ) and

f δ (x, u) := f (x, u) if x ∈ Ω 1 δ (u -ϕ(x)) if x ∈ B \ Ω
and as before w = det(D 2 u) -1 and U is the cofactor matrix of D 2 u. In view of (4.2) and the definition of f δ , it is natural to introduce the functional defined over convex functions on B by

J δ 1 (v) := Ω F (x, v(x))dx + 1 2δ B\Ω (v -ϕ) 2 -F B (v)
where

F B (v) := B log(det(∂ 2 v))
so that

J δ 1 (v) = J 1 (v) + 1 2δ B\Ω (v -ϕ) 2 - B\Ω log(det(∂ 2 v)). (4.3)

A priori estimates for the second boundary value problem

0 ≥ (F ′ B (u) -F ′ B ( u))(u -u) = B (U ij w ij -U ij w ij )(u -ϕ) + ∂B ψ(U ij -U ij )∂ i (u -ϕ)ν j = Ω f (x, u)(u -ϕ) - B f (u -ϕ) + 1 δ B\Ω (u -ϕ) 2 + ∂B ψ(U νν -U νν )∂ ν (u -ϕ) ≥ Ω f (x, ϕ)(u -ϕ) - B f (u -ϕ) + 1 δ B\Ω (u -ϕ) 2 + ∂B ψ(U νν -U νν )∂ ν (u -ϕ)
where, in the last line, we have used the fact that ∇u -∇ϕ = ∂ ν (uϕ)ν on ∂B and set

U νν = Uν • ν, U νν = U ν • ν.
Using the fact that f (x, ϕ), f , ϕ, ∇ϕ and U are bounded, we thus get

1 δ B\Ω (u -ϕ) 2 + ∂B ψU νν ∂ ν u ≤ C 1 + B |u| + ∂B ∂ ν u + ∂B U νν . (4.5)
Denoting by R the radius of B and by the same argument as in Lemma 2.2 in [START_REF] Chau | Monge-Ampère functionals and the second boundary value problem[END_REF], one has

U νν = 1 R d-1 ∂ ν u d-1 + E with |E| ≤ C(1 + ∂ ν u d-2 ) on ∂B. (4.6)
Moreover since u is convex and u = ϕ = 0 on ∂B, one has max

B |u| = -min B u ≤ 2R∂ ν u(x) for all x ∈ ∂B. (4.7)
Putting together (4.5), (4.6), (4.7) and the fact that inf ∂B ψ > 0, we obtain

∂B (∂ ν u) d ≤ C(1 + ∂B (∂ ν u) d-1 )
which gives a bound on ∂ ν u L d (∂B) hence also on max B |u| by (4.7) so that finally the bound on δ -1 B\Ω (uϕ) 2 follows from the latter bounds and (4.5).

Existence and uniqueness of a smooth uniformly convex solution

Thanks to Theorem 1.1 in [START_REF] Chau | Monge-Ampère functionals and the second boundary value problem[END_REF], a Leray-Schauder degree argument and the a priori estimate (4.4), one easily deduces the following:

Theorem 4.2. For every δ > 0, the second boundary value problem (4.2) admits a unique uniformly convex solution which is W 4,p (B) for every p ∈ [1, +∞).

Proof.

Let D := {u ∈ C(B), u C(B) ≤ C + 1}
where C is the constant from (4.4). For t ∈ [0, 1] and u ∈ D, it follows from Theorem 1.1 in [START_REF] Chau | Monge-Ampère functionals and the second boundary value problem[END_REF] that there exists a unique W 4,p for every p ∈ [1, ∞) and uniformly convex solution of

V ij w ij = tf δ (x, u), w = det(D 2 v) -1 in B, v = ϕ, w = ψ on ∂B (4.8)
where V denotes the cofactor matrix of D 2 v. We denote by v = T t (u) the solution of (4.8). Moreover, by Theorem 2.1 of [START_REF] Chau | Monge-Ampère functionals and the second boundary value problem[END_REF], for every α ∈ (0, 1) there are a priori bounds on v C 3,α and on sup B (det Finally, uniqueness follows from the same argument as in Lemma 7.1 from [START_REF] Trudinger | Boundary regularity for the Monge-Ampère and affine maximal surface equations[END_REF] where it is proven that two smooth solutions actually have the same gradient on ∂B and then are the minimizers of the same strictly convex minimization problem hence coincide.

(D 2 v) + det(D 2 v) -1 ) that only depend on C, α, δ, ϕ C 3,1 and the constant λ in (2.1). Therefore (t, u) ∈ [0, 1] × D → T t (u) is continuous on [0, 1] × D and T t is compact in C(B)
In dimension d = 2, following the argument of Remark 4.2 of Trudinger and Wang [START_REF] Trudinger | The affine Plateau problem[END_REF] and taking advantage of the fact that the right-hand side of the Abreu equation (4.2) does not depend on δ on Ω, we have the following local bound (which we have used in the proof of Proposition 3.4): We also have

z ij = w ij w - w i w j w 2 + 2 δ ij η + 2 η i η j η 2 -u ijk u k -u ik u jk , (4.10)
multiplying by wU = [D 2 u] -1 , using

wU ij z ij ≥ 0 at x 0 , U ij w ij = f (x, u) ≤ C
and the identities

wU ij u ik u jk = u ii = ∆u, wU ij u ijk u k = - w k w u k = - ∇w w • ∇u, (4.11)
(the second identity is classically obtained by first differentiating the relation -log(w) = log(det D 2 u) and then taking the scalar product with ∇u) as well as the fact that Tr(U) = ∆u in dimension d = 2, we get

0 ≤ C -wU ∇w w • ∇w w + 2 w η ∆u + 2wU ∇η η • ∇η η + ∇w w • ∇u -∆u. (4.12)
Using (4.9) and using again wU = [D 2 u] -1 , we then obtain

wU ∇w w • ∇w w = 4wU ∇η η • ∇η η + D 2 u∇u • ∇u + 4 ∇η η • ∇u (4.13) and ∇w w • ∇u = 2 ∇η η • ∇u + D 2 u∇u • ∇u. (4.14)
Replacing (4.13), (4.14) in (4.12), multiplying by η and rearraging gives (4.15) gives η(x 0 )∆u(x 0 ) ≤ 2C ′ and since ∆u(x 0 )w(x 0 ) 1/2 ≥ 2 we get the desired lower bound on the minimum of η -2 w. In the remaining case w(x 0 ) ≥ 1 4 η(x 0 ) ≥ η 2 (x 0 ) 2r 2 and we reach the same conclusion. This gives a local lower bound on w i.e. the desired local upper bound on det(D 2 u).

∆u(η -2w) ≤ Cη -2 w η U∇η •∇η -2∇η •∇u ≤ Cη + ∇η L ∞ ∇u L ∞ (Ω) ≤ C ′ (4.15) If η(x 0 ) ≥ 4w(x 0 ),

Convergence

Letting δ → 0 for fixed ε

In this paragraph, we fix ε (and thus normalize it to ε = 1 as we did in the whole of section 4).

Theorem 5.1. Let u δ be the unique smooth strictly convex solution of (4.2), then u δ converges uniformly on Ω to the unique minimizer of (4.1) as δ → 0 + . Proof. We already know from (4.4) that (possibly up to an extraction) u δ converges locally uniformly on B to some convex u and it also follows from (4.4) that u ∈ S[ϕ, Ω]. Let v ∈ S[ϕ, Ω] (extended by ϕ on B \ Ω), thanks to (4.2) and the convexity of J δ 1 we first have We now have to pay attention to the boundary term, we know from (4.6) that θ δ := ψU νν δ satisfies 0 ≤ θ δ ≤ C(1 + (∂ ν u δ ) d-1 ) so that thanks to (4.4), θ δ is bounded in L d d-1 (∂B), up to an extraction we may therefore assume that it weakly converges in L d d-1 (∂B) to some nonnegative function θ. By convexity we also have that for τ > 0

J δ 1 (v) -J δ 1 (u δ ) ≥ ∂B U νν δ ψ∂ ν (u δ -ϕ) i.e. J 1 (v) -J 1 (u δ ) ≥ 1 2δ B\Ω (u δ -ϕ) 2 + B\Ω (log(det(D 2 ϕ)) -log(det(D 2 u δ ))) + ∂B U νν δ ψ∂ ν (u δ -ϕ)
∂ ν u δ (x) ≥ D τ,ν u δ (x) := 1 τ u δ (x -τ ν(x)) -u δ (x -2τ ν(x)) , ∀x ∈ ∂B
For small fixed τ > 0 note that D τ,ν u δ is bounded independently of δ thanks to (4.4) and that it converges as δ → 0 pointwise to D τ,ν ϕ, we thus have lim inf

δ→0 ∂B θ δ ∂ ν (u δ -ϕ) ≥ lim inf δ→0 ∂B θ δ (D τ,ν u δ -∂ ν ϕ) = ∂B θ(D τ,ν ϕ -∂ ν ϕ)
where in the last line we have passed to the limit using the fact that we have the product of a weakly convergent sequence with a strongly convergent sequence. Now letting τ → 0 and using the smoothness of ϕ, we deduce that lim inf

δ→0 ∂B U νν δ ψ∂ ν (u δ -ϕ) ≥ 0.
Since J 1 is lower semi-continuous thanks to Lemma 2.2, we can conclude that

J 1 (v) ≥ lim inf δ→0 J 1 (u δ ) ≥ J 1 (u)
so that u solves (4.1) and by uniqueness of the minimizer there is in fact convergence of the whole sequence.

In the previous proof we have used:

Lemma 5.2. Let u δ be the unique smooth strictly convex solution of (4.2) as before, then

lim sup δ→0 B\Ω log(det(D 2 u δ )) ≤ B\Ω log(det(D 2 ϕ)).
Proof. The key point here is the estimate B ∆u δ = ∂B ∂ ν u δ ≤ C which follows from (4.4). Let ω be an arbitrary Borel subset of B, we have (for some constant C varying from a line to another):

ω log(det(D 2 u δ )) ≤ C(|ω| + ω det(D 2 u δ ) 1/2d ) ≤ C |ω| + ω ∆u δ ≤ C |ω| + |ω| 1/2 B ∆u δ 1/2 = C |ω| + |ω| 1/2 ∂B ∂ ν u δ 1/2 so that ω log(det(D 2 u δ )) ≤ C(|ω| + |ω| 1/2 ). ( 5.1)
Take 0 < R ′ < R with Ω contained in B R ′ (recall R is the radius of B), we then have, thanks to Lemma 2.2, the fact that log(det(D 2 ϕ)) is bounded and (5.1):

lim sup δ→0 F B\Ω (u δ ) = lim sup δ→0 F B\Ω (u δ ) ≤ lim sup δ→0 F B R ′ \Ω (u δ ) + lim sup δ→0 F B\B R ′ (u δ ) ≤ F B R ′ \Ω (ϕ) + C(|B \ B R ′ | + |B \ B R ′ | 1/2 ) ≤ F B\Ω (ϕ) + C ′ (|B \ B R ′ | + |B \ B R ′ | 1/2 ).
The desired result follows by letting R ′ tend to R.

Full convergence

We now take δ = δ ε > 0 with

lim ε→0 + δ ε = 0, ( 5.2)
i.e. we only have a single small parameter ε and we consider the secondboundary value problem

εU ij ε w ε ij = g ε (x, u ε ), in B, u ε = ϕ, w ε = ψ on ∂B (5.3)
where ψ := det((D 2 ϕ) -1 ),

g ε (x, u) := f (x, u) if x ∈ Ω 1 δε (u -ϕ(x)) if x ∈ B \ Ω , w ε = det(D 2 u ε ) -1
and U ε is the cofactor matrix of D 2 u ε . We further assume that there is an α > 0 such that

(f (x, u) -f (x, v))(u -v) ≥ α(u -v) 2 , ∀(u, v) ∈ R d ,
and a.e. x ∈ Ω (5.4) which amounts to say that the integrand F is uniformly convex in its second argument. Under these assumptions, we have a full convergence result:

Theorem 5.3. Let u ε be the unique smooth strictly convex solution of (5.3), then u ε converges uniformly on Ω to the unique minimizer of (1.1) as ε → 0 + .

Proof.

Step 1: a priori estimates. The first step of the proof is similar to the proof of Proposition 4.1. Again define u := ϕ, U as the cofactor matrix of D 2 ϕ, w := det(D 2 ϕ) -1 and f ε := ε U ij w ij . We then have together with (5.4):

0 ≥ ε(F ′ B (u ε ) -F ′ B ( u))(u ε -u) ≥ Ω (f (x, ϕ) -f ε )(u ε -ϕ) + α Ω (u ε -ϕ) 2 + 1 δ ε B\Ω (u ε -ϕ) 2 + ε ∂B ψ(U νν ε -U νν )∂ ν (u ε -ϕ)
thanks to the fact that f (x, ϕ)f ε is bounded uniformly with respect to ε, using Young's inequality and invoking (4.6), we get

Ω (u ε -ϕ) 2 + 1 δ ε B\Ω (u ε -ϕ) 2 + ε ∂B (∂ ν u ε ) d ≤ C.
(5.5)

Step 2: convergence. Thanks to (5.5), up to taking a subsequence of vanishing ε n , we may assume that u ε converges locally uniformly in B to some u such that u = ϕ in B \ Ω so that the restriction of u to Ω belongs to S[ϕ, Ω]. For every v convex on B such that v = ϕ on ∂B, define Let then v ∈ S[ϕ, Ω] (extended by ϕ on B \ Ω), we then have

J ε (v) -J ε (u ε ) ≥ ε ∂B ψU νν ε ∂ ν (u ε -ϕ) hence J 0 (v) ≥ lim inf ε J 0 (u ε ) + lim inf ε ε(F B (v) -F B (u ε )) -lim sup ε ε ∂B ψU νν ε ∂ ν ϕ.
Arguing as in the proof of Proposition 3.5, we may actually assume that F Ω (v) > -∞ so that lim inf ε εF B (v) ≥ 0. As for an upper bound for εF B (u ε ) we use the fact that thanks to (5.5), we have ∂B ∂ ν u ε ≤ Cε

Lemma 2 . 2 .

 22 If ω is an open subset of B with ω ⊂⊂ B, then for every sequence of convex functions u n converging locally uniformly on B to u, one has lim sup n

Following a similar convexity 4 . 1 .

 41 argument as in Lemma 2.2 in Chau and Weinkove[START_REF] Chau | Monge-Ampère functionals and the second boundary value problem[END_REF], we first have Proposition Let u be a smooth and uniformly convex solution of (4.2), then maxB |u| + ∂B |∂ ν u| d + 1 δ B\Ω |u -ϕ| 2 ≤ C (4.4)for some constant C only depending on B, ϕ C 3,1 (B) and the constant λ in (2.1).Proof. First observe that by convexity and (2.1), u < 0 in B and ∂ ν u > 0 on ∂B. Define u := ϕ, U as the cofactor matrix of D 2 ϕ, w := det(D 2 ϕ) -1 and f := U ij w ij (whose L ∞ norm only depends on ϕ C 3,1 (B) and the constant λ in (2.1)) we have by the concavity of F B , (4.2) and the monotonicity of f (x, .):

  for every t ∈ [0, 1]. Since T 0 is constant and by(4.4) it has a unique fixed point in D, again by (4.4), T t has no fixed point on ∂D, it thus follows from the Leray-Schauder Theorem that T 1 has a fixed point in D, this proves the existence claim for (4.2).

Proposition 4 . 3 .

 43 Let d = 2 and u be the solution of (4.2) then for every open set ω ⊂⊂ Ω, det(D 2 u) L ∞ loc (ω) is bounded independently of δ.Proof. Let B r := B r (0) ⊂⊂ Ω, and observe that thanks to (4.4) both f δ (., u(.)) L ∞ (Ω) = f (., u(.)) L ∞ (Ω) and ∇u L ∞ (Ω) are bounded independently of δ. Define then η(x) := 1 2 (r 2 -|x| 2 ) and consider z := log(w) -2 log(η) -1 2 |∇u| 2 , by construction z achieves its minimum at an interior point x 0 of B r at such a point, we have

  det(D 2 ϕ)) -log(det(D 2 u δ ))) + ∂B U νν δ ψ∂ ν (u δϕ). It follows from Lemma 5.2 below that lim inf δ→0 B\Ω (log(det(D 2 ϕ)) -log(det(D 2 u δ ))) ≥ 0.

JF

  ε (v) := Ω (x, v(x))dx + 1 2δ ε B\Ω (vϕ) 2ε B log(det(∂ 2 v)).

  -1/d and argue in a similar way as in the proof of Lemma 5.2, to obtainεF B (u ε ) ≤ Cε(1 + B det(D 2 u ε ) 1/d ) ≤ Cε(1 + ∂B ∂ ν u ε ) ≤ C(ε + ε 1-1/d ), which yields lim inf ε ε(F B (v) -F B (u ε )) ≥ 0.14 Thanks to (4.6), we have∂B ψU νν ε ∂ ν ϕ ≤ C ∂B (1 + (∂ ν u ε ) d-1 )but, thanks to (5.5) and Hölder's inequality, we deduceε ∂B (∂ ν u ε ) d-1 ≤ Cε 1 d so that J 0 (v) ≥ lim inf ε J 0 (u ε ) = J 0 (u)hence u solves (1.1) (and the whole family u ε converges uniformly on Ω to u by uniqueness of the minimizer of J 0 on S[ϕ, Ω]).

Acknowledgments: G.C. is grateful to the Agence Nationale de la Recherche for its support through the project MAGA (ANR-16-CE40-0014).


Aybüke Öztürk
email: aybuke.ozturk@univ-lyon2.fr

Stéphane Lallich
email: stephane.lallich@univ-lyon2.fr

Jérôme Darmont
email: jerome.darmont@univ-lyon2.fr

A Visual Quality Index for Fuzzy C-Means

Keywords: Fuzzy Clustering, Fuzzy C-Means, Quality Indices, Visual Index, Elbow Rule

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Clustering refers to the assignment of unlabeled data points into clusters (groups) so that the points belonging to the same cluster are more similar to each other than those within different clusters. There are various types of clustering strategies, including crisp and fuzzy clustering. In crisp (or hard) clustering, a data point can belong to one and only one cluster, while in fuzzy clustering [START_REF] Ruspini | Numerical methods for fuzzy clustering[END_REF], a data point can belong to several clusters. Fuzzy clustering is very useful in many applications, e.g., the text categorization of various news into different clusters: a science, a business, and a sport cluster; where an article containing the keyword "gold" could belong to all three clusters. Furthermore, it is also possible to open discussions with domain experts when using fuzzy clustering.

Clustering algorithms behave differently for different reasons. The first reason relates to dataset features such as geometry and the density distribution of clusters. The second reason is the choice of input parameters such as the fuzziness coefficient m (m = 1 indicating that clustering is crisp and m > 1 that clustering becomes fuzzy).

These parameters all affect the quality of clustering. To study how the choice of parameters impacts clustering quality, we need a quality criterion. For instance, when the dataset is well separated and has only two variables, a scatter plot can help determine the number of clusters. However, when the dataset has more than two variables, a good quality index is needed to compare various cluster configurations and choose the appropriate number of clusters.

Achieving a good clustering involves both minimizing intra-cluster distance (compactness) and maximizing inter-cluster distance (separability). A common issue in this process is that clusters are split up while they could be more compact. Many cluster quality indices have been proposed to address this problem for hard and fuzzy clustering, but none of them is always highly efficient [START_REF] Pal | Correction to" on cluster validity for the fuzzy c-means model[END_REF].

Moreover, there is no real-life golden standard for clustering analysis, since various experts may have different points of views about the same data and express different constraints on the number and size of clusters. Thanks to a visual index, different solutions can be presented with respect to the data. Thus, experts can make a trade-off between their opinion and the best local solutions proposed by the visual index.

Hence, in this paper, we first review existing quality indices that are wellsuited to fuzzy clustering, such as [START_REF] Bezdek | Cluster validity with fuzzy sets[END_REF][START_REF] Chen | Rule-base self-generation and simplification for datadriven fuzzy models[END_REF][START_REF] Caliński | A dendrite method for cluster analysis[END_REF][START_REF] Fukuyama | A new method of choosing the number of clusters for the fuzzy c-mean method[END_REF][START_REF] Xie | A validity measure for fuzzy clustering[END_REF][START_REF] Zhang | A novel cluster validity index for fuzzy clustering based on bipartite modularity[END_REF]. Then, we propose an innovative, visual quality index for the well-known Fuzzy C-Means (FCM) method. Moreover, we compare our proposal with state-of-the-art quality indices from the literature on several numerical real-world and artificial datasets.

The remainder of this paper is organized as follows. Section 2 recalls the principles of fuzzy clustering. Section 3 surveys quality indices for fuzzy clustering. Section 4 details our visual quality index. Section 5 reports on the experimental comparison of our quality index against existing ones on different datasets. Finally, we conclude this paper and provide research perspectives in Section 6.

Principles of Fuzzy Clustering

Fuzzy inertia is a core measure in fuzzy clustering. Fuzzy inertia F I (Equation 1) is composed of fuzzy within-inertia F W (Equation 2) and fuzzy between-inertia F B (Equation 3). Membership coefficients u ik of data point i to cluster k are usually stored in a membership matrix U that is used to calculate F W , F B and F I. Note that F I = F W + F B. Moreover, F I is not constant because it depends on u ik . When F W changes, the values of F I and F B also change.

F I = n i=1 K k=1 u m ik d 2 (x i , x) (1)
F W = n i=1 K k=1 u m ik d 2 (x i , c k ) (2)
F B = n i=1 K k=1 u m ik d 2 (c k , x) (3)
where n is the number of instances, K is the number of clusters, m is the fuzziness coefficient (by default, m = 2), c k is the center of the k th cluster ∀1 ≤ k ≤ K, x is the grand mean (the arithmetic mean of all data -Equation 4), and function d 2 () computes the squared Euclidean distance.

x = 1 n n i=1 x i (4)
FCM is a common method for fuzzy clustering that adapts the principle of the K-Means algorithm [START_REF] Macqueen | Some methods for classification and analysis of multivariate observations[END_REF]. FCM, proposed by [START_REF] Dunn | A fuzzy relative of the isodata process and its use in detecting compact well-separated clusters[END_REF] and extended by [START_REF] Bezdek | Fcm: The fuzzy c-means clustering algorithm[END_REF], applies on numerical data. Since numerical data are the most common case, we choose to experiment our proposals with FCM.

The aim of the FCM algorithm is to minimize F W . It starts by choosing K data points as initial centroids of the clusters. Then, membership matrix values u ik (Equation 5) are assigned to each data point in the dataset. Centroids of clusters c k are updated based on Equation 6 until a termination criterion is reached successfully. In FCM, this criterion can be a fixed number of iterations t, e.g., t = 100. Alternatively, a threshold can be used, e.g., = 0.0001. Then, the algorithm stops when

F W K+1 / | F W K+1 -F W K | < . u ik = 1 K j=1 ( xi-c k 2 xi-cj 2 ) 1 m-1 (5) c k = n i=1 u m ik x i n i=1 u m ik (6)
3 Fuzzy Clustering Quality Indices

According to Wang et al. [START_REF] Wang | On fuzzy cluster validity indices[END_REF], there are two groups of quality indices. Quality indices in the first group are based only on membership values. They notably include partition coefficient index V P C [3] (Equation 7; 1 K ≤ V P C ≤ 1; to be maximized) and Chen and Linkens' index V CL [START_REF] Chen | Rule-base self-generation and simplification for datadriven fuzzy models[END_REF] (Equation 8; 0 ≤ V CL ≤ 1; to be maximized). V CL takes into consideration both compactness (first term of V CL ) and separability (second term of V CL ).

V P C = 1 n n i=1 K k=1 u 2 ik (7) V CL = 1 n n i=1 max k (u ik ) - 1 c K-1 k=1 K j=k+1 1 n n i=1 min(u ik , u ij ) , (8)
where

c = K-1 k=1 k.
Quality indices in the second group associate membership values to cluster centers and data. They include an adaptation of the Ratio index V F Ratio to fuzzy clustering [START_REF] Caliński | A dendrite method for cluster analysis[END_REF] (Equation 9; 0 ≤ V F Ratio ≤ +∞; to be maximized), Fukuyama and Sugeno's index V F S [START_REF] Fukuyama | A new method of choosing the number of clusters for the fuzzy c-mean method[END_REF] (Equation 10; -F I ≤ V F S ≤ F I; to be minimized), and Xie and Beni's index V XB [START_REF] Xie | A validity measure for fuzzy clustering[END_REF][START_REF] Pal | On cluster validity for the fuzzy c-means model[END_REF] (Equation 11; 0 ≤ V XB ≤ F I/n * min x j -v k 2 ; to be minimized).

V F Ratio = F B/F W (9) V F S = F W -F B (10)
V XB = K k=1 n i=1 u m ik x i -v k 2 n * min j,k v j -v k 2 (11)
When the number of clusters increases, the value of quality indices mechanically increases, too. Then, the important question is: how useful is the addition of a new cluster? To answer this question, the most common solutions are penalization and the Elbow Rule [START_REF] Cattell | The scree test for the number of factors[END_REF].

The first way to penalize a quality index is to multiply it by a quantity that diminishes the index when the number of clusters increases. In this case, the main difficulty is to choose the penalty. For instance, the penalized version of V F Ratio is Calinski's V F CH [START_REF] Caliński | A dendrite method for cluster analysis[END_REF] (Equation 12; 0 ≤ V F CH ≤ +∞; to be maximized), where the penalty is based on both the number of clusters and data points.

V F CH = F B/(K -1) F W/(n -K) = n -K K -1 F B F W (12)
The second way to penalize a quality index is to evaluate index evolution relatively to the number of clusters, by considering the curve of the index' successive values. The most appropriate value of K can be determined visually by help of the Elbow Rule or algebraic calculation [START_REF] Dimitriadou | An examination of indexes for determining the number of clusters in binary data sets[END_REF].

To construct a visual determination of the Elbow Rule, K is represented on the horizontal axis and the considered quality index on the vertical axis. Then, we look for the value of K where there is a change in the curve's concavity. This change represents the optimal number of clusters K. To construct an algebraic determination, let i K being the index value for K clusters. The variation of i K before K and after K are compared. In case of a positive Elbow, the second difference min K ((i

K+1 -i K ) -(i K -i K-1
)) is minimized. Yet, since the values before K and after K are used for calculation, the Elbow Rule can be applied to more than two clusters only. Among all the above-stated quality indices, there is no single quality index that gives the best result for any dataset. Thus, there is room for a new quality index that is specifically tailored for fuzzy validation and helps the user choose the value of K.

An Index Associated with a Visual Solution

Building a new quality index, we first consider F W to evaluate compactness and F B to evaluate separability. We can choose to calculate either F B -F W , which is similar to V F S except for the sign, or F B ÷ F W , which is similar to V Adding a new cluster often improves clustering quality mechanically. Thus, many authors penalize the quality index with respect to K (the smaller n is, the greater the penalty), e.g., V F CH (Section 3). To obtain a penalized index, SF D is first linearly transformed in an index belonging to [0, 1], obtaining the Transformed Standardized Fuzzy Difference T SF D (Equation 13; T SF D ∈ [0, 1]; to be maximized). Finally, by penalizing T SF D as V F CH , we obtain the Penalized Standardized Fuzzy Difference P SF D (Equation 14; P SF D ∈ [0, (n -K)/(K -1)]; to be maximized).

T SF D = 1 + SF D 2 = F B F I (13)
P SF D = T SF D * n -K K -1 = F B -F W F I * n -K K -1 (14)
Instead of penalizing the quality index, another solution is to visualize the search for the best number of clusters K. First solution is to apply the Elbow Rule to T SF D. T SF D is plotted with respect to K in Figure 1(a). The drawback of this method is that the horizontal axis corresponds to an arithmetic scale of K values, which is not satisfying. To fix this problem, we suggest to plot F B with respect to F I, which we call Visual T SF D. Our aim is not to give an automatic solution, but to help the user visually choose the most appropriate K value. The visualization we propose is shown in Figure 1(b), where the blue line plots T SF D with respect to K, the full red line is the diagonal that corresponds to the best solutions (F B = F I) such that T SF D = 1, and the dashed red line connects the origin to each point associated with K values. The smaller the angle between the full red line and the dashed red line, the better is the solution. As the value of K increases, the angle between the dashed red line and the diagonal decreases. Then, we choose the value of K beyond which the decrease becomes negligible. This value is considered as the optimal number of clusters. For example, in Figure 1(b), a first solution could be K = 4, a better solution K = 6, and it is not very interesting to consider K > 6.

Experimental Validation

In this section, we compare our proposals T SF D, P SF D, Visual TSFD and the use of the Elbow Rule to state-of-the-art clustering quality indices for FCM-like clustering algorithms, i.e., V P C , V CL , V F CH , V F S and V XB (Section 3).

In our experiments, the FCM algorithm is parameterized with its default settings: termination criterion = 0.0001 and default fuzziness coefficient m = 2. All clustering quality indices are coded in Python version 2.7.4. 1)

Datasets

Quality indices are compared on ten real-life datasets (Table 1; IDs 1-10) from the UCI Machine Learning Repository1 and seven artificial datasets (Table 1; . In real-life datasets, the true number of clusters is assimilated to the number of labels. Although using the number of labels as the number of clusters is debatable, it is acceptable if the set of descriptive variables explains the labels well. In artificial datasets, the number of clusters is known by construction. Moreover, we created new artificial datasets by introducing overlapping and noise to some of the existing datasets, such as E1071-3 [START_REF] Meyer | Package e1071[END_REF], Ruspini [START_REF] Ruspini | Numerical methods for fuzzy clustering[END_REF] and E1071-5 [START_REF] Meyer | Package e1071[END_REF] (Table 1; IDs 12-14). To create a new dataset, new data points are introduced, and each must be labeled. To obtain a dataset with overlapping, we modify the construction of the E1071 artificial datasets [START_REF] Meyer | Package e1071[END_REF]. In the original datasets, there are three or five clusters of equal size (50). Cluster i is generated according to a Gaussian distribution N (i; 0.3). To increase overlapping in the three clusters while retaining the same cluster size, we only change the standard deviation from 0.3 to 0.4. Then, there is no labeling problem. To introduce noise in a dataset, we add in each cluster noisy points generated by a Gaussian variable around each label gravity center. Noisy data are often generated by distributions with positive skewness. For example, in a two-dimensional dataset, for each label, we add points that are far away from the corresponding gravity center, especially on the right-hand side, which generally contains the most points. Then, we draw a random number r between 0 and 1. If r ≤ 0.25, the point is attributed to the left-hand side. Otherwise, the point is attributed to the right-hand side. This method helps obtain noisy data that are 1 / 4 times smaller and 3 / 4 times greater, respectively, than the expected value for the considered label. This process is applied to the Ruspini dataset [START_REF] Ruspini | Numerical methods for fuzzy clustering[END_REF].

Experimental Results

In our experiments, all validation indices (Sections 3 and 4) are applied on all the datasets from Table 1. Moreover, since presenting all the results would take too much space, we retain only the best results for each index (even excluding P SF D). As shown in Table 1, it is more difficult to predict an appropriate number of clusters for real-life datasets than for artificial datasets. Considering all indices, the average rate of success is indeed 21% in the case of real data, against 66% in the case of artificial data. Whatever the type of data, Visual TSFD outperforms the other indices, with 5 wins out of 10 in the case of real datasets, and 6 wins out of 7 in the case of artificial datasets. The worst results are obtained with V P C and V XB (0/10 and 4/7 wins each). The other indices achieve intermediary results. In addition, when the value given by Visual TSFD is erroneous, it is quite close to the expected K, in contrast to V F S , our closest competitor (Table 1; Wine, Glass, Segmentation, Ecoli and Bensaid). For example, the optimal number of clusters should be 6 for the Glass dataset. V F S = 12, Visual TSFD's results are 5 and 7. Furthermore, we compare in Figures 2 and3 Visual TSFD and the plot obtained with the Elbow Rule (which is labeled Elbow T SF D) with respect to K, on a sample of both real-life and artificial datasets bearing different characteristics, i.e., Glass, Vehicle, Ecoli, Ruspini, Ruspini noised and E1071-5overlapped (Table 1). As is clearly visible from Figures 2 and3 value. This makes our method more flexible that the existing ones in real-life scenarios.

Conclusion and Perspectives

In this paper, we propose a novel quality index for FCM called Visual TSFD, which provides an overview of fuzzy clustering with respect to the number of clusters. We compare Visual TSFD to several clustering quality methods from the literature and experimentally show that it outperforms existing methods on various datasets. Furthermore, Visual TSFD can also be used in the case of categorical data with Fuzzy K-Medoids [START_REF] Park | A simple and fast algorithm for k-medoids clustering[END_REF]. Thus, Visual TSFD allows to deal with heterogeneous datasets, which makes our method a simple but noteworthy contribution, in our opinion. As a result, our next step is to design an ensem-ble fuzzy clustering method based on Visual TSFD that would deal with both numerical and categorical data.

  F Ratio . Unfortunately, F I = F B +F W is not constant and F B -F W ∈ [-F I, +F I]. To take this particularity of fuzzy clustering into account, we propose to standardize F B -F W by considering the Standardized Fuzzy Difference SF D = (F B -F W ) ÷ F I instead. Then, SF D ∈ [-1, +1].

Fig. 1 :

 1 Fig. 1: Comparison of Elbow Rule (a) and Visual TSFD (b) on the WineQuality-Red dataset (Table 1)

  , Visual TSFD gives a better visual idea than Elbow T SF D. Elbow T SF D indeed highlights K values of 3 or 4, while the T F SD blue plot systematically indicates larger K values.

Fig. 2 :Fig. 3 :

 23 Fig. 2: Comparison of Elbow T SF D and Visual TSFD (1/2)

Table 1 :

 1 Quality Indices Experiment Results with Different Datasets

	ID	Datasets	# of data points	# of clusters	VP C VCL F B VF CH VF S VXB	Elbow VT SF D	Visual VT SF D
	1	Wine	178	3	2	2	8	12	8	2	3	5
	2	Iris	150	3	2	2 3	3	3	2	3	3
	3	Seeds	210	3	2	3 3	3	3	2	3	3
	4	Glass	214	6	2	2 12 12	12 2	4	5,7
	5	Vehicle	846	4	2	2	2	2	5	2	3	4,5
	6	Segmentation	2310	7	2	4	4	4	12 12	3	7,8
	7 Movement Libras 360	15	2 18 16 16	18 2	14	14,16
	8	Ecoli	336	8	2	3	3	3	12 3	3	3,7
	9	Yeast	1484	10	2	2	5	2	12 2	4	7,8
	10 WineQuality-Red 1599	6	2	2 6	7	6	2	3	6
	11	Bensaid [17]	49	3	3	3 9	11	11 3	3	5
	12	E1071-3 [16]	150	3	3	3 3	3	3	3	3	3
	13	Ruspini [1]	75	4	4	4 4	4	4	4	3	4
	14	E1071-5 [16]	250	5	2	5 4	5	5	2	3	5
	15	E1071-3-overlapped	150	3	2	3 3	2	3	2	3	3
	16 Ruspini noised	95	4	4 12 4	4	4	4	4	4
	17	E1071-5-overlapped	250	5	2	2	4	5	4	2	3	5
	# of wins for real-life datasets 0	1	3	2	3	0	3	5
	# of wins for artificial datasets 4	5	4	5	5	4	4	6
		Total # of wins		4	6	7	7	8	4	7	11

http://archive.ics.uci.edu/ml/

Acknowledgments

This project is supported by the Rhône Alpes Region's ARC 5: "Cultures, Sciences, Sociétés et Médiations" through A. Öztürk's Ph.D. grant.