Ripon K Chakrabortty
email: ripon.chakrabortty@student.adfa.edu.au

Ruhul A Sarker

Daryl L Essam

Event Based Approaches for Solving Multi-mode Resource Constraints Project Scheduling Problem

Keywords: Multi-mode Resource Constrained Project Scheduling Problem, Mixed Integer Linear Programming formulations, Event based formulations

Over the last few decades, a number of mathematical models have been introduced for solving Multi-mode Resource Constrained Project Scheduling Problems (MRCPSPs). However the computational effort required in solving those models depends on the number of variables. In this paper, we attempt to reduce the number of variables required in representing MRCPSPs by formulating two new event-based models. A comparative study was conducted by solving standard benchmark instances using a common objective function for the developed as well as the existing mathematical models. The study provided interesting insights about the problem characteristics, model sizes, solution quality, and computational effort of these approaches.

Introduction

The Resource Constrained Project Scheduling Problem (RCPSP) has been a challenging research topic over the last few decades [START_REF] Hartmann | A survey of variants and extensions of the resourceconstrained project scheduling problem[END_REF]. For it, each activity in RCPSP requires some resources for its processing, such as: machines and tools, human and their skills, raw materials, natural resources (energy, water, land, etc.), information and money. A variant of RCPSP is known as the Multi-mode Resource Constrained Project Scheduling Problem (MRCPSP) and for it each activity can be performed using one out of a set of different resources (known as modes), with specific activity durations and resource requirements. MRCPSP provides the most comprehensive framework for realistically representing the flexibilities that decision makers have [START_REF] Zapata | The multimode resource constrained multiproject scheduling problem: Alternative formulations[END_REF].

MRCPSPs have many applications within computer integrated manufacturing, planning in make-to-order production, as well as other production planning problems, such as flow shop and job shop scheduling problem. These applications have stimulated many mathematical formulations for the MRCPSPs. From the last decades or so, CPM and PERT were the fundamental techniques to solve project scheduling problems that minimize the project completion time, while satisfying the precedence rela-tionships among the activities, with an assumption that the resources are unlimited. In reality, as the resources are limited, the formulation of MRCPSPs becomes complex and the formulation varies with the assumptions made.

To solve RCPSPs, integers linear programming (ILP) based approaches are widely used [START_REF] Koné | Comparison of mixed integer linear programming models for the resource-constrained project scheduling problem with consumption and production of resources[END_REF]. In the formulation by Pritsker et al. [START_REF] Pritsker | Multiproject scheduling with limited resources: A zero-one programming approach[END_REF], the binary decision variable is defined to be one if activity finishes at time instant , and to be 0 otherwise. This formulation requires the use of at most nT binary decision variables and ( + ) constraints, where n is the number of jobs and m is the number of arcs. Later, Klein and Scholl [START_REF] Klein | Computing lower bounds by destructive improvement: An application to resource-constrained project scheduling[END_REF] developed two different models where binary variable is defined to be one if activity is in progress in period , and to be 0 otherwise. These formulations require nT binary decision variables and (

) constraints, which is more than the previous one. Therefore, real life projects with many activities are impossible to solve by using exact algorithms, as they involve a large number of variables and constraints. Again, according to Lawler [START_REF] Lawler | Combinatorial optimization: networks and matroids[END_REF], the MRCPSP belongs to ; and the solution time is a function of ( ), which is further increased exponentially with the modes, number of variables and constraints [START_REF] Sabzehparvar | A mathematical model for the multi-mode resource-constrained project scheduling problem with mode dependent time lags[END_REF]. Many researchers have tested their approaches by solving benchmark instances from PSPLIB (Project Scheduling Problem Library), with different number of activities, such as 30, 60, 90 and 120 for single mode and 10-30 activities for multi-modes. According to Horroelen and Leus [START_REF] Herroelen | Project scheduling under uncertainty: Survey and research potentials[END_REF] and Alcaraz et al. [START_REF] Alcaraz | Solving the Multi-Mode Resource-Constrained Project Scheduling Problem with genetic algorithms[END_REF], many of the 60-activity and most of the 90-and 120-activity instances from PSPLIB are beyond the capability of exact methods. Koné et al. [START_REF] Koné | Event-based MILP models for resource-constrained project scheduling problems[END_REF] has recently indicated that MILP model can deal with up to 25-35 activities within the convenient computational time of 500 seconds. For this reason, any mathematical model that can reduce the number of variables and constraints in representing MRCPSPs would be very practically useful.

Pinto and Grossmann [START_REF] Pinto | A continuous time MILP model for short term scheduling of batch plants with pre-ordering constraints[END_REF] provided an example from a plastic compound plant where the activity duration varied from 0.783 to 11.25 days. The conventional discrete-time formulations are not suitable to tackle this class of problem unless the processing time is rounded or the time horizon is decomposed into smaller time intervals, which only yields approximate solutions. Although Koné et al. [START_REF] Koné | Event-based MILP models for resource-constrained project scheduling problems[END_REF] presented two event based approaches for single mode RCPSP that can deal with non-integer duration, there is no such model available for MRCPSP. In this regard, the continuous time formulation, along with the event based formulation, can be interesting approaches to study MRCPSPs. From the literature review, it is clear that the trend in MRCPSP research and practice is to formulate the problem as a mathematical model, such as Mixed Integer Linear Programming (MILP), and then solve the model using a conventional optimization technique. However, considering the state of current computational capability, the conventional methods are incapable of solving MILP models with a large number of variables and constraints. In some models, the number of variables increases exponentially with the number of time periods, which encourages the development of alternative mathematical models which by redefining the variables, attempt to minimize the number of variables and constraints.

In this paper, we have formulated two event based models for representing MRCPSPs. For reducing the number of binary decision variables, in those models, we have considered continuous time or time-independent approaches that are known as 'Start/End Event based formulation (SEE) and 'On/Off Event based formulation (OOE)'. The proposed formulations are variants of that of Zapata et al. [START_REF] Zapata | The multimode resource constrained multiproject scheduling problem: Alternative formulations[END_REF], but with fewer variables. Thus, instead of the 3 binary variables for Zapata et al. [START_REF] Zapata | The multimode resource constrained multiproject scheduling problem: Alternative formulations[END_REF], it took 2 variables for SEE and binary variables for OEE, where m represents number of modes and n represents the number of activities. It was conjectured from the earlier researchers that the solution time of integer programming can be bounded as an exponential function of the number of variables [START_REF] Lenstra | Integer programming with a fixed number of variables[END_REF]. So, the main intention of this paper is to reduce the number of variables and also the solution times for solving MRCPSP. Our developed mixed integer programming models were solved by using a commercial optimization solver. The proposed model is capable of dealing with MRCPSPs of up to 30-35 activities within reasonable computing time (500 sec). To demonstrate the applicability of the developed models, a number of 10-30 activity multi-mode benchmark instances from the Project Scheduling Library (PSLIB) were solved. The solutions were analyzed for problems with different levels of complexity and are compared with the traditional discrete time approach. The computational studies were also conducted to analyze the effects of different factors that relate to performance.

The structure of the paper is as follows: in section 2, we define the basic MRCPSPs. The terminologies, the proposed event based MILP formulations are described in section 3. In section 4, solution approaches are discussed. The experimental studies, along with the computational results and analysis are provided in section 5. Finally, we provide conclusions in the last section.

2

Problem Description MRCPSP belongs to the set of combinatorial optimization problems [START_REF] Koné | Event-based MILP models for resource-constrained project scheduling problems[END_REF]. The MRCPSP under study in this paper is based on the following assumptions: (i) the activities composing a project have certain and known durations that pertain to certain modes; all predecessors must be finished before an activity can start; (iii) resources can be either renewable or both renewable and non-renewable. Renewable resources are available in limited amounts, whereas non-renewable resources are fixed for the whole project; (iv) activities are non-preemptive (i.e., cannot be interrupted when in progress); (v) the main objective is to minimize the project completion time. Let I be the number of activities to be scheduled, R be the number of available renewable resources to be allocated and be the number of modes available for activity i. The activities constituting the project are represented by a set {0,..., I+1}, where 0 and I+1 are the dummy nodes representing the start and end respectively. Accordingly, the set of renewable resources are defined by {0,1,...,r}, whereas the non-renewable resources are defined by {0,1,...,w}. The other important notations are depicted in the next section as Nomenclature. The traditional formulations for scheduling are discrete in manner, i.e., it depends on discretization of the time horizon, T, which is represented as {0, ..., t}. Each resource type has a certain capacity limit which cannot be exceeded throughout the project life.

The Mathematical Models

In this section, firstly, we have presented a discrete time approach and then two new formulations are developed. The standard MRCPSP requires sequencing the project activities, so that the precedence constraints are met, the execution mode for each activity is determined, the resource constraints are met and the project duration is minimized. In spite of having several mathematical models for MRCPSP [START_REF] Zapata | The multimode resource constrained multiproject scheduling problem: Alternative formulations[END_REF] and [START_REF] Beşikçi | Different resource management policies in multimode resource constrained multi-project scheduling[END_REF], the model from Talbot [START_REF] Talbot | Resource-constrained project scheduling with time-resource tradeoffs: The nonpreemptive case[END_REF] is still being used as a basis for discrete time approach and hence is employed throughout this paper as the discrete time (DT) approach for comparing with our proposed models. In the MRCPSP model, a binary decision variable is defined to be 1 if activity starts at mode m at time instant , and 0 otherwise. The model can be presented as follows:

Minimize Z = ∑ ∑ ∈ ∈ Constraints: ∑ ∑ = 1 ∀ ∈ (1) ∑ ∑ tx ≥ ∑ ∑ (t + d )x (i, j) ∈ P (2) ∑ ∑ ∑ ≤ ∀ ∈ ∀ ∈ ∈ (3) ∑ ∑ ∑ ≤ ∀ ∈ ∀ ∈ ∈ (4) = ∑ ∑ ( + ) ∈ ∈ ∀ ∈ (5) ∈ {0,1} ∀ ∈ , ∀ ∈ , = 1, … , (6)
Here, constraint equation [START_REF] Hartmann | A survey of variants and extensions of the resourceconstrained project scheduling problem[END_REF] represents that every job or activity must be handled exactly one time for all modes. Constraint set (2) ensures that precedence relationships are maintained. Meanwhile constraint sets (3) and ( 4) respectively represent that the capacity of the renewable and non renewable resources are satisfied. Finally, constraint set (5) defines the finish time of activity i. , is used for the quantity of resource k that is required immediately after event e at mode m. That is why this proposed approach is termed as a start/end event based formulation (SEE).

Min

= 0 (7)

≥ + - 1 - ∀( , ) ∈ , > , ∀ ∈ , ∀ ∈ (8)
≥ ∀ ∈ , < (9)
∑ ∑ = 1 ∀ ∈ (10)
∑ ∑ = 1 ∀ ∈ (11)
∑ ∑ + ∑ ∑ ≤ 1 ∀( , ) ∈ , ∀ ∈ (12) = ∑ ∑ ∈ ∀ ∈ (13) = , + ∑ ∑ -∑ ∑ ∈ ∈ ∀ ∈ , ≥ 1, ∈ (14)
≤ ∀ ∈ , ∈ , ∀ ∈ (15)
≤ ≤ LS x + LS (1 -x ) ∀i ∈ A, ∀e ∈ ε, ∀m ∈ M (16) ≤ ≤ (17) ( + ) ≤ ≤ ( + ) + (1 - ) ∀ ∈ , ∀ ∈ , ∀ ∈ (18) ≥ 0 ∀ ∈ (19) ∑ ∑ ∑ x ≤ ∀ ∈ ∀ ∈ ∈ (20) ≥ 0 ∀ ∈ , ∈ , ∀ ∈ (21) ∈ {0,1}, ∈ {0,1} ∀ ∈ ∪ {0, + 1}, ∀ ∈ , ∀ ∈ (22)
Constraint [START_REF] Sabzehparvar | A mathematical model for the multi-mode resource-constrained project scheduling problem with mode dependent time lags[END_REF] stipulates that event 0 starts at time 0. Inequalities [START_REF] Herroelen | Project scheduling under uncertainty: Survey and research potentials[END_REF] ensure that if activity starts at event e and ends at event f, then ≥ + . Any other combination of values for and yield either ≥ or ≥ -, which are redundant with constraint [START_REF] Alcaraz | Solving the Multi-Mode Resource-Constrained Project Scheduling Problem with genetic algorithms[END_REF]. Constraints [START_REF] Koné | Event-based MILP models for resource-constrained project scheduling problems[END_REF] and [START_REF] Pinto | A continuous time MILP model for short term scheduling of batch plants with pre-ordering constraints[END_REF] ensure that the start and end events must have a single occurrence for all of its available modes respectively. Constraint [START_REF] Lenstra | Integer programming with a fixed number of variables[END_REF] describes the precedence relationship between activities. Constraint [START_REF] Beşikçi | Different resource management policies in multimode resource constrained multi-project scheduling[END_REF] represents that the total resource demands for each resource starts at event 0. Constraint [START_REF] Talbot | Resource-constrained project scheduling with time-resource tradeoffs: The nonpreemptive case[END_REF] defines that for each resource r, its demand immediately after event e is equal to its demand immediately after the previous event e-1, plus the demand required by the activities (on any mode m) that start at event e, minus the demand required by the activities that end at event e. Constraint (15) limits the demanded resources at each event to the availability of resources. Constraints ( 16)-( 18) are valid inequalities based on activity time windows. Inequality (20) is for non-renewable resource w, while the usages of that resource throughout the whole event date must to be within its maximum available limit.

3.4

On/Off Event Based Continuous Time Formulation (OOE)

In this section, we propose another model that uses only one type of binary variable per event. Similar to the last section, this is also an extension of an earlier model proposed by Koné et al. [START_REF] Koné | Event-based MILP models for resource-constrained project scheduling problems[END_REF]. The SEE model needs more variables for representing events than doe the OEE model [START_REF] Koné | Comparison of mixed integer linear programming models for the resource-constrained project scheduling problem with consumption and production of resources[END_REF]. For it, the number of events is exactly equal to the number of activities. The representation of the resource constraints is also simpler and easier. Again, rather than using tree decision variables for start/end formulation, here in this model only two decision variables are used. A decision variable is set to be 1 if activity starts at event e at mode m or if it is still being processed immediately after event e. Here, in the same way as in the previous model, a continuous decision variable ( ) represents the date of event e. The OOE formulation is presented below.

Min

∑ ∑ ∈ ≥ 1 ∀ ∈ (23) ≥ + - , ∀ ∈ , ∀ ∈ , ∀ ∈ (24) = 0 (25) ≥ ∀ ≠ -1 ∈ (26) ≥ + - , - - , -1 ∀( , , ) ∈ × , > , ∀ ∈ (27) ∑ ∑ ≤ 1 - - , ∀ ∈ \{0} (28) ∑ ∑ ≤ ( -) 1 + - , ∀ ∈ \{0} (29) + ∑ ∑ ≤ 1 + (1 - ) ∀ ∈ , ∀( , ) ∈ (30) ∑ ∑ ≤ ∀ ∈ , ∀ ∈ (31) ≤ ≤ - , + 1 - - , ∀ ∈ , ∀ ∈ , ∀ ∈ (32) ≤ ≤ ∀ ∈ , ∀ ∈ (33) ≥ 0 ∀ ∈ (34) ∈ {0,1} ∀ ∈ , ∀ ∈ , ∀ ∈ (35)
The objective is once again to minimize the makespan. Constraints (23) ensure that each activity is processed at least once during the project. Constraints (24) build up the relationship between the makespan and the event dates: ≥ + if is in process at event e but not at event e-1, i.e., if starts at event e. Constraints (25) and (26) ensure the basic definition of event sequencing. Constraints (27) link the binary optimization variable to the continuous variable , and ensure that if activity starts immediately after event e and ends at event f, then the date of event f is at least equal to the date of event e plus the processing time of activity . Constraints (28) and (29) are called contiguity constraints, and they ensure non pre-emption (i.e. that the events after which a given activity is being processed are adjacent). Constraints (30) maintain the precedence relationship. Constraints (31) are the resource constraints limiting the total demand of activities in process at each event. Constraints (32) and (33) ensure that the start time of any activity is between its earliest and latest start time.

Experimental Study

The mathematical models developed in the previous section for MRCPSPs are integer programming models which can be solved using standard optimization algorithms. The developed models were coded with the LINGO optimization software, and were executed on an Intel core i7 processor with 16.00 GB RAM and a 3.40 GHz CPU. For solving the models, the Branch and Bound algorithm (B&B) algorithm within LINGO was applied. To compare the models, we have selected a number of MRCPSP bench-mark instances from the popular test library known as PSPLIB-Project Scheduling Problem Library [START_REF] Kolisch | Characterization and generation of a general class of resource-constrained project scheduling problems[END_REF]. The characteristics of the MRCPSP instances vary a lot in terms of the average number of resources required (resource factor, RF), the resource feature incorporates time parameters (resource strength, RS) and finally the integration of resource and precedence features. For justifying the proposed formulations, we considered the RS and RF for both renewable and non renewable resource to be 0.5, which means that if the project has 4 renewable or non renewable resources then for any particular time each activity can demand two resources.

As suggested by the earlier models, the conventional DT approach needs binary variables, where m is the number of modes, n is the number of activities including the dummy ones and T is the total project completion time. But for the SEE and OOE approaches, they need only 2 and binary variables to represent any MRCPSP. This indicates the binary variables required for SEE are twice that of OOE. To compare the above mentioned three models, we have considered a number of 10activity test instances, with low activity duration time, that have an average makespan of 18. As shown in Table 1, the conventional DT approach needs on average 648 binary variables, while OOE and SEE require 432 and 864 respectively. After multiplying each activity duration time by 10 (now considered as a higher activity duration), the number of binary variables for DT becomes 6156 and the makespan increases from 18 to 81. However, the number of variables for SEE and OOE remain the same. As the number of variables is the dominating factor for the computational time required in solving a model, the time recorded in Table 1 consistently shows that the models with higher numbers of variables takes longer computational time than the same ones with lower numbers of variables. We have further analyzed and compared the above mentioned three models using a number of test instances from PSLIB with 10, 12, 14, 16, 18, 20 and 30 activities for both high and low duration times. Here, high duration means five times longer than low duration. For the low activity duration, the number of variables and the computational time required, with respect to the number of activities, are presented in Figures 1 and2 respectively. For the high activity duration, these plots are shown in Figures 3 and4. In these figures, for either low or high activity duration of all three models, the number of variables increases with the number of activities and the computational time consistently increases with the increase of the number of variables. There is no impact on the number of variables in SEE and OOE when the activity duration is changed (say changing from low to high). However, as can be seen in Table 1 and Figures 1 and3, the number of variables significantly increases in DT as activity du-ration increases. From these four figures, it is clear that OOE is the best model irrespective of the length of activity duration. However, SEE is better than DT only for high activity duration. In real life projects, some activities may be added or removed during project execution. In this section, we will analyze the impact of new activity insertion. To do this, we assumed that a set of new activities with precedence relations need to be added to the project network. Here we have considered three cases for a 10-activity instance (J1010_4) with randomly generated duration and resource usage requirements on each mode. In figure 5, Ins-2 represents the insertion of two new activities which replaced 2 old activities. From the figure, it is quite clear that the impact of new activity insertion is insignificant for event based approaches. However it has a clear impact on DT model. Note that, for all approaches, the number of constraints increases significantly with the increase of new activities. In a project, the activity priority may be revised and this may change the precedence relationship. In this section, we have also analyzed the impacts of any particular precedence modifications. An example of precedence modification is when the project network for A needs to satisfy the additional precedence relationship , and that it also no longer needs to satisfy the relations in . We have considered four randomly generated precedence modifications for analyzing the impact on the three models. As shown in figure 6, similar conclusion can be drawn for precedence modifications as was earlier found for new activity insertions. From the figure, it can be observed that with increasing numbers of precedence modifications (PM), the number of variables remain same for both event based approaches, while for the DT approach, the number of variables is increased. So, it can be claimed, that the proposed event based approaches are independent of precedence modifications. But for all approaches, the number of constraints is slightly increased with the increase of precedence modifications.

Fig. 6. Distribution of number of variables with precedence modifications [START_REF] Klein | Computing lower bounds by destructive improvement: An application to resource-constrained project scheduling[END_REF]

Conclusions

In this paper, we have developed two new event based mathematical models (mixed integer linear program) for MRCPSP that reduced the number of variables and hence the required solution times. For solving the models, the branch and bound techniques were implemented within LINGO. A number of benchmark problems with different numbers of activities and resources were solved and compared. From the comparisons, it can be concluded that (i) there is no impact on the number of variables in SEE and OOE with activity duration, (ii) the number of variables significantly increases in DT with increases of activity duration, (iii) OOE is the best model irrespective of the length of activity duration, (iv) SEE is better than DT only for high activity durations, (v) the effect of new activity insertion or precedence modification on the number of variables are insignificant for SEE and OOE, and (vi) there is a clear negative effect for new activity insertion or precedence modification on the DT model. So practitioners will benefit by choosing the proposed event based model (OOE) for MRCPSP, as it requires the lower number of variables and lower computational time in comparison to other existing models. The future research should be on proposing some real time event based RCPSP approaches for handling uncertainties or research disruptions. Implementing this model for more complex problems with large number of activities could be another important extension of these proposed models.

Fig. 1 .Fig. 2 .Fig. 3 .Fig. 4 .

 1234 Fig. 1. Relationship between number of variables and activities with low duration time

Fig. 5 .

 5 Fig. 5. Distribution of number of variables with new activity insertion

  In discrete time formulations, the variables are indexed by time. In contrast to that formulation, in this section, a new event based formulation is proposed for MRCPSP. It is based on the concept of the single mode RCPSP model proposed by Koné et al.[START_REF] Koné | Event-based MILP models for resource-constrained project scheduling problems[END_REF]. For classical MRCPSP, the numbers of activities are often considered as n+2 (where n is the activity number with 0 & n+1 being the dummy source and sink nodes). On the contrary, this event based formulation considers only (n+1) activities as events. Event-based formulations also have the advantage that they can deal with non-integer activity processing times. More importantly, for instances with long scheduling horizons, event-based models involve fewer variables in comparison to the models indexed by time. Here, we used only two types of binary variables one for start event () at mode m and the other one is for the end event ( ) at mode m. A continuous variable represents the date of event e and a continuous ble,

	3.3	Start/End Event Based Continuous Time Formulations (SEE)

Table 1 .

 1 Binary variable and constraint based comparison among all mentioned approaches

	10-activity Instances	Parameters	DT	SEE	OOE
	With low activity-	No. of Variables	648	864	432
	duration time	Solution time (sec)	2.0	3.0	1.5
	With high activity-	No. of Variables	6156	864	432
	duration time	Solution time (sec)	5.0	3.0	1.5


Juan Carlos Figueroa-García

Emanuel Medina
email: ejmedinap@gmail.com

David Rubio-Espinosa

Non-cooperative Games involving Type-2 Fuzzy Uncertainty: An approach

Classical games theory is based on a set of deterministic payoffs that represent different strategies taken by its players. Deterministic (crisp) payoffs are used by homogeneous players while in many social and business scenarios, different uncertainties such as numerical, linguistic, stochastic, etc are involved. This way, we proposed a model for a two-player game problem involving nonprobabilistic uncertainty coming from uncertain fuzzy payoffs, through Interval Type-2 fuzzy numbers. A method for solving this kind of uncertain games is proposed, and an introductory example is shown.

Introduction and Motivation

Game theory became an important discipline, widely applied in fields such as economics, financial, engineering, etc (see Nash [START_REF] Nash | Non cooperative games[END_REF], and Neumann and Morgenstern [START_REF] Neumann | Theory of Games and Economic Behavior[END_REF]). In some cases, uncertainty affects the payoffs of a game, so in real scenarios both the behavior of the payoffs and the strategies taken by players become an important topic in game theory. In classical (crisp) two-person games the payoffs are deterministic numbers, but in some cases the payoffs contain imprecision that can be handled through fuzzy numbers; in this case the game is called a Fuzzy Game.

Fuzzy games in literature are based on the use of classical fuzzy sets (a.k.a Type-1 fuzzy sets) alongside Linear Programming (LP) methods. The most important works have been proposed by Bector and Chandras [START_REF] Bector | Fuzzy Mathematical Programming and Fuzzy Matrix Games[END_REF], Campos [START_REF] Campos | Fuzzy linear programming models to solve fuzzy matrix games[END_REF], Delgado, Verdegay and Vila [START_REF] Delgado | Playing Matrix Games Defined by Linguistic Labels[END_REF], Li [START_REF] Li | An effective methodology for solving matrix games with fuzzy payoffs[END_REF], Butnariu [START_REF] Butnariu | Fuzzy games a description of the concept[END_REF][START_REF] Butnariu | Solution concept for n-Person games[END_REF], Vijay et al [START_REF] Vijay | Bimatrix games with fuzzy payoffs and fuzzy goals[END_REF], Monroy et al [START_REF] Monroy | Set-valued cooperative games with fuzzy payoffs. the fuzzy assignment game[END_REF], Cunlin and Qiang [START_REF] Cunlin | Nash equilibrium strategy for fuzzy non-cooperative games[END_REF], and Larbani [START_REF] Larbani | Non cooperative fuzzy games in normal form: A survey[END_REF]. All those works handle imprecision around the payoffs in the game using Type-1 fuzzy sets/numbers.

It is commonly assumed that the payoffs of a game are perfectly defined without uncertainty. In some cases the payoffs cannot be surely defined due to multiple uncertainties, and they are defined by multiple people who express their perceptions about the payoffs of a game. As usual, those perceptions are not deterministic values and they involve multiple aspects that every person believes possible, so different people's thoughts lead to have multiple fuzzy sets. This way, a fuzzy game composed by multiple people defining its payoffs needs special tools to handle this kind of fuzzy uncertainty.

Our proposal is based on the use of Type-2 fuzzy numbers to involve all perceptions and opinions of who define the payoffs of a game, dealing with non-consensus and linguistic uncertainty around their perceptions. Using LP methods, the extension principle, and the decomposition theorem, we provide a mathematical framework for solving a matrix game conformed by two-players under uncertainty coming from disagreement and non-consensuses on the definition of the payoffs of the game.

We also use interval Type-2 fuzzy numbers to represent the set of possible solutions that the players can use, given a particular perception about the payoffs matrix. A method for solving such games, using α-cuts is proposed. The method is explained through an example.

This paper is divided into five sections. Section 1 shows the Introduction and Motivation. In Section 2, the classical two-players game is described, and its LP formulation is presented. In Section 3, some basics on Type-2 fuzzy numbers are recalled. In Section 4, a Type-2 fuzzy game model is presented. Section 5 presents the proposal for solving the problem; Section 6 presents an application example, and finally in Section 7 some concluding remarks are suggested.

The classical two-players game

In a two-players game, they are trying to make a move to earn the best possible benefits, given possible moves of the other player and known payoffs of every move. In zero-sum games, a move gives a gain to a player while its opponent losses the same. Then a game G is conformed by a set S m of strategies of the Player I and a set S n of strategies of the Player II. If Player II selects the pure strategy s j and the Player I selects the pure strategy s i , then Player II has to pay a ij units to Player I.

All payoffs of the game are comprised into a payoff matrix namely A, where each element a ij is the payoff of the Player I, and a negative value means that Player I has to pay a ij to Player II, as shown next:

A =           a 11 a 12 • • • a 1j • • • a 1n a 21 a 22 • • • a 2j • • • a 2n . . . . . . . . . . . . . . . . . . a i1 a i2 • • • a ij • • • a in . . . . . . . . . . . . . . . . . . a m1 a m2 • • • a mj • • • a mn          
where

a ij ∈ R, i ∈ {1, 2, • • • , m} and j ∈ {1, 2, • • • , n}.
This matrix A was originally defined with deterministic (crisp) payoffs, which leads to have only one strategy that maximizes the game. Decision making for Player I uses a maxmin principle which seeks to maximize its income through the worse possible movings done by its opponent, so the final decision consists in doing a move i that earns the maximum of the minimum payoffs given all possible moves j of its opponent (maxmin principle is applied to Player II).

LP formulation

Basically, an LP formulation of a two-players game is intended to find the optimal frequency of moves which maximizes the value of the game. To do so, Player I has to select a move x * i (S m ) to obtain the best expected payoff (a ij ).

Max x∈S m Min m i=1 a i1 x i , m i=1 a i2 x i , • • • , m i=1 a in x i s.t. m i=1 x i = 1 (1)
x i 0

Analogous model can be formulated for Player II. Dantzig [START_REF] Dantzig | Linear Programming and Extensions[END_REF] proposed the first LP model for games based on the works of Nash [START_REF] Nash | Non cooperative games[END_REF], and Neumann & Morgenstern [START_REF] Neumann | Theory of Games and Economic Behavior[END_REF]. His proposal uses an auxiliary variable v that represents the min operator. The LP formulation for Player I is shown next.

Max z = v s.t. v - m i=1 a ij x i 0 ∀ j ∈ S n (2)
m i=1 x i = 1 x i 0, v ∈ R v is an auxiliary variable which operates as v = Min ( m i=1 a i1 x i , • • • , m i=1 a in x i ).
The LP model for Player II is the dual model of Player I, so the optimal value of the game is the same for both players. Further details can be found in Dantzig [START_REF] Dantzig | Linear Programming and Extensions[END_REF].

Type-2 fuzzy sets

A Type-2 fuzzy set (see Mendel [START_REF] Mendel | Uncertain Rule-Based Fuzzy Logic Systems: Introduction and New Directions[END_REF]) is an ordered pair {(x, µ Ã(x)) : x ∈ X}, where A is its linguistic label and à represents uncertainty around the word A. It is defined as follows:

à = {((x, u), µ Ã(x, u) | ∀x ∈ X, ∀u ∈ J x ⊆ [0, 1]}
where J x is the set of primary memberships of Ã, 1 µ Ã(x, u) 0, and u is its domain of uncertainty.

An Interval Type-2 fuzzy set (IT2FS) (see Figure 1) is a special kind of T2FS in which µ Ã(x, u) = 1. Hence, an IT2FS Ã is composed by an infinite amount of Type-1 fuzzy sets embedded into its primary memberships J x . In other words (see Mendel, John & Liu [START_REF] Mendel | Interval type-2 fuzzy logic systems made simple[END_REF]) Here, Ã is an Interval Type-2 fuzzy set defined over a domain (usually R + ), its support supp( Ã) is enclosed into the interval a ∈ [A L(+) , A R(+) ]. µ ã has a linear Type-2 fuzzy set with parameters A L(+) , A L(-) , A R(-) , A R(+) and Ā. FOU is the Footprint of Uncertainty of the Type-2 fuzzy set and A e is an embedded Type-1 fuzzy set into its FOU. Uncertainty A is defined by its FOU as the union of all primary memberships, i.e.

à = {(x, µ Ã(x)) | x ∈ X} where µ Ã(x) is completely characterized by J x ⊆ [0, 1]. A Ae UMF LMF 1 µ à A ∈ R A L(+) A L(-) A R(-) A R(+) Ā F OU ( Ã)
FOU( Ã) = x∈X J x (3)
The FOU of à involves all J x weighted by a secondary membership function f x (u)/(x, u), in this case 1/(x, u). Any FOU is bounded by two membership functions: an Upper membership function (UMF)= µ à and a Lower membership function (LMF)= µ à which have e embedded sets (A e ).

The Zadeh's extension principle

The selected approach to to model fuzzy functions is the Zadeh's Extension principle (see Klir and Yuan [START_REF] Klir | Fuzzy Sets, Uncertainty and Information[END_REF] which is presented as follows Theorem 31 Let f be a function such as f :

X 1 , X 2 , • • • , X n → z, and A i is a fuzzy set in X i , i = 1, 2, • • • , n with x i ∈ X i , then we have f (A 1 , A 2 , • • • , A n )(z) = sup z=f (x1,x2,••• ,xn) min i {µ A1 (x 1 ), µ A2 (x 2 ), • • • , µ An (x n )} (4)
Then, we can project any function z = f (x 1 , x 2 , • • • , x n ) into a fuzzy set using their membership functions µ A1(x1) , µ A2(x2) , • • • , µ An(xn) . For simplicity, we compute the set µ z using a mapping of A i (x i ) instead of x i . This is by using α-cuts and the representation theorem of fuzzy set (See Klir and Yuan [START_REF] Klir | Fuzzy Sets and Fuzzy Logic: Theory and Applications[END_REF]), presented next.

The extended representation theorem

The α-cut of A is defined as

α A = {x | µ A (x)
α}. The function z = f (x * ) comes from f ( Ãi )(z), so f ( Ãi )(z) can be easily computed through α-cuts instead of mapping x i . Then, we compute a function of f -1 (µ Ã) → x to finally compute functions of fuzzy sets. The extension of the α-cut of A to the α-cut of à has been proposed by Figueroa [START_REF] Figueroa-García | An approximation method for type reduction of an interval Type-2 fuzzy set based on α-cuts[END_REF] as the union of all possible α-cuts over A e , α A e = {x | µ Ae (x) α}:

α Ã = {(x, u) | J x α} ∀ u ∈ [0, 1]
Now, the boundaries of each α-cut are defined as follows (see Figure 2):

α Ã L ij = α Ã L(+) ij , α Ã L(-) ij (5)
α Ã R ij = α Ã R(-) ij , α Ã R(+) ij (6) 1 µ
à A ∈ X μ à µ à 1 à α à α ÃR(-)α ÃR(+) α ÃL(-) α ÃL(+) F OU

Fig. 2. Fuzzy set Ã

Using (4), we can extend z to a Type-2 fuzzy environment, as follows:

f ( Ã)(z) = f α ( Ã1 , Ã2 , • • • , Ãn )(z) = f ( α Ã1 (x 1 ), α Ã2 (x 2 ), • • • , α Ãn (x n ))(z) (7)

Type-2 fuzzy numbers (IT2FN)

A Type-2 fuzzy number (T2FN) can be seen as the extension of a Type-1 fuzzy number. Then à is a T2FS whose UMF and LMF are normalized fuzzy numbers, α A is a closed interval for all α ∈ [0, 1], and bounded support supp(A). Now, if both the UMF and LMF of à are fuzzy numbers, then à is an IT2FN. Figure 1 depicts a Type-2 fuzzy number since both its UMF and LMF are fuzzy numbers as well.

In some applications, a crisp measure (seen either as the expected or most possible value) of an IT2FS, is required. The computation of the expected value of an IT2FN is possible using a centroid-based method, whose details can be found in Figueroa [START_REF] Figueroa-García | An approximation method for type reduction of an interval Type-2 fuzzy set based on α-cuts[END_REF], Karnik & Mendel [START_REF] Karnik | Centroid of a type-2 fuzzy set[END_REF], and Melgarejo [START_REF] Duran | Improved iterative algorithm for computing the generalized centroid of an interval type-2 fuzzy set[END_REF].

The IT2FS games model

We have based our proposal on the Li's α-cuts method for fuzzy games (see Li [START_REF] Li | An effective methodology for solving matrix games with fuzzy payoffs[END_REF]). In our paper, the payoffs matrix has been defined using Interval Type-2 fuzzy uncertainty, so the problem (in terms of Player I) is as follows:

Max x∈S m Min m i=1 ãi1 x i , m i=1 ãi2 x i , • • • , m i=1 ãin x i (8)
s.t. m i=1 x i = 1 (9)
x i 0

In this model each ãij is defined as an IT2FN, so its mathematical programming form (for Player I), is:

Max z = ṽ s.t. ṽ - m i=1 ãij x i 0 ∀ j ∈ S n (11) m i=1
x i = 1

x i 0

Every α Ã leads to four optimal strategies per cut. So every boundary of α Ã reaches an optimal value namely α ZL(+) , α Z L(-) , α ZR(-) and α ZR(+) :

α Ã L(+) → α Z L(+) (12) α Ã L(-) → α Z L(-) (13) α Ã R(-) → α Z R(-) (14) α Ã R(+) → α Z R(+) (15)
After computing all values of the game regarding α then the fuzzy set of optimal solutions namely Z is composed by using f -1 ( α Ã)(z).

Solving IT2FGs

The problem of solving IT2FLPs has been investigated by Figueroa [START_REF] Figueroa-García | An approximation method for type reduction of an interval Type-2 fuzzy set based on α-cuts[END_REF][START_REF] Figueroa | A general model for linear programming with interval type-2 fuzzy technological coefficients[END_REF][START_REF] Figueroa | Computing optimal solutions of a linear programming problem with interval type-2 fuzzy constraints[END_REF][START_REF] Figueroa-García | A transportation model with interval type-2 fuzzy demands and supplies[END_REF], so the optimization of IT2FGs become a natural extension of his results. Note that we face a problem of having infinite Type-1 fuzzy sets into the FOU of Ã, so we have to deal with the FOU of ãij instead to solve infinite fuzzy LPs. Now, the proposed method to compute a solution of IT2FG modeled as an LP when having a Type-2 fuzzy payoff matrix is described as follows:

Proposed α-cut based method 1. Select k ∈ N amount of α-cuts 2.
Compute all the α-cuts of à 3. Solve the 4k LP related problems regarding ( 12), ( 13), ( 14) and ( 15) 4. Compose f ( Ã)(z) using equation ( 7) 5. Compose Z using f -1 ( α Ã)(z) 6. Compute C( Z), z l and z u using Figueroa [START_REF] Figueroa-García | An approximation method for type reduction of an interval Type-2 fuzzy set based on α-cuts[END_REF] 7. Return [z l , z u ] as the interval of expected values of the game 8. Defuzzify Z using c z = (z l + z u )/2

Decision Making

Fuzzy Decision making is not based on having a single (global) optimum of the problem, they are based on the idea of having a set of optimal solutions coming from fuzzy information, so our problem needs to compute a set of possible solutions of the game called Z, as a function of the initial fuzzy payoffs ãij .

Every particular value of α ãij leads to a particular value of the game z * which is satisfactory in different degrees J z ∈ [µ Z , μ Z ], so the idea is to cover as much as possible scenarios. The centroid of Z provides an idea of what an expectation of the game is, so the analyst can use it to have an idea of the expected value of the game.

Finally, what it is proposed is a map of all possible values of a ij in order to provide the analyst a better idea of the behavior of the fuzzy game.

Application example

The example is based on a classical example introduced by Lieberman [START_REF] Hillier | Introduction to Operations Research[END_REF]. The original example is a crisp 4 × 5 two-players zero-sum game, where the solution is easily computed using an LP model. Now suppose there is no consensus about the payoffs of the game, so different uncertainties affect its definition (i.e ambiguity, multiple people's opinions, imprecision, etc).

The analyst only has crisp payoffs, but they are not reliable. This way, some experts of the game are asked for their opinions about its payoffs, so we have to comprise all the information given by the experts using Type-2 fuzzy numbers which collects linguistic uncertainty (coming from experts opinions and perceptions) around a ij to finally provide a solution. The initial crisp payoffs are shown next: 1. Crisp payoffs of the example For simplicity, experts were asked to express their opinion around the crisp payoffs shown in Table 1 using the statement "What is your opinion about optimistic and pessimistic payoffs of the strategy (i, j)?". The words Optimistic and Pessimistic are used to compose triangular fuzzy numbers (namely T (a, b, c)) for every UMF and LMF. Note that a means pessimistic, b means crisp (see Table 1), and c means optimistic payoff. The collected IT2FNs are presented in Table 2.

(i,j) j = 1 j = 2 j = 3 j = 4 j = 5 i = 1 1 -3 2 -2 1 i = 2 2 3 0 3 -2 i = 3 0 4 -1 -3 -2 i = 4 4 0 -2 2 -1 Table
Lower membership functions (UMF) of the example (i,j) T(-6,-4,2) T(-1,0,3) T(-6,-2,2) T(-2,2,5) T(-4,-1,1.5) Table 2. Parameters of the example Our example uses the values shown in Table 2 to compose the matrix à whose elements are µ ãij . Now, the game is solved using the α-cuts decomposition and the representation theorem.

j = 1 j = 2 j = 3 j = 4 j = 5 i = 1 T(0.7,1,1.

Computational results

We have selected k = 21 equally distributed α-cuts and computed α Ã using ( 5) and [START_REF] Li | An effective methodology for solving matrix games with fuzzy payoffs[END_REF]. After solving the 4k LP problems related to (1) and (2) (this is, an LP for (12), ( 13), [START_REF] Mendel | Uncertain Rule-Based Fuzzy Logic Systems: Introduction and New Directions[END_REF], and ( 15)), we compose f ( Ã)(z) using step 3) and equation [START_REF] Butnariu | Fuzzy games a description of the concept[END_REF]. The set of optimal solutions Z is shown in Figure 3. A clearer clear idea of the behavior of the game is given by the centroid of Z, so we compute z l and z u using Figueroa [START_REF] Figueroa-García | An approximation method for type reduction of an interval Type-2 fuzzy set based on α-cuts[END_REF], whosse results are shown next:

C( Z) ∈ [z l , z u ] = [-0.17637, 0.47788] → c z = 0.15075
Now, C( Z) is expected interval of values that both groups of players can expect. Different strategies can be taken in the game, for which we provide a set of possible solutions x * i and y * j (Player I and II respectively) that fits into C( Z). This is, a set of solutions whose values z * fit into C( Z) having different membership degrees. Note that as more α-cuts are computed, we have more strategies into C( Z).

Discussion of the results

We have obtained different solutions coming from different payoffs. For instance, the value z * = 0.427 is within C( Z) and comes from 0.8 ÃL(-) . The optimal strategy (probabilities) for Player I is x * = [0.554, 0.425, 0.021, 0], and the optimal strategy for Player II is y * = [0, 0.274, 0, 0.097, 0.629]. Other example comes from 1 ÃR(+) which is equivalent to the crisp solution of the game (outside C( Z)); its optimal value is z * = -0.306, the optimal strategy for Player I is x * = [0.565, 0.354, 0.081, 0], and the optimal strategy for Player II is y * = [0, 0.291, 0, 0.048, 0.661].

Other choice is to select a particular payoff. Suppose that players select a payoff, for instance 0.8 ÃL(-) . This selection leads to z * = 0.427 whose set of memberships is μ Z = 0.8 and µ Z ≈ 0.9057. Note that the value of μ Z is the solution of the game for 0.8 ÃR(-) , and the value µ Z is computed using the Type-reduction algorithm proposed by Figueroa [START_REF] Figueroa-García | An approximation method for type reduction of an interval Type-2 fuzzy set based on α-cuts[END_REF], so the solution satisfies both groups of players in the interval α ∈ [ 0.8, 0.9057 ].

Concluding Remarks

We have presented a method to solve games involving Type-2 fuzzy uncertainty, using the results presented by Li [START_REF] Li | An effective methodology for solving matrix games with fuzzy payoffs[END_REF] combined with IT2FNs. It is possible to deal with this kind of uncertainty using well known algorithms, so its applicability to real scenarios is promissory.

The α-cuts approach is a map of the optimal solutions of the IT2FG, so as more complete the map of à is, the more complete Z is, giving more information to the analyst. Although there is no a single solution of the IT2FG, we provide a map of the possible strategies that players can take in different scenarios, which is important in decision making under uncertainty.

The use of natural language in optimization and decision making has opened a door to use hybrid methodologies such as the presented in this paper. Different approaches can be proposed to solve this kind of games, so our proposal is just a method for modeling and solving the game.

Fig. 1 .

 1 Fig. 1. Interval Type-2 Fuzzy set Ã

1Fig. 3 .

 3 Fig. 3. Set of optimal values of the game Z

⋆ Juan Carlos Figueroa-García is Assistant Professor at the Universidad Distrital Francisco José de Caldas, Bogotá -Colombia. ⋆⋆ Emanuel Medina is undergraduate student of the Industrial Engineering Dept. of the Universidad Distrital Francisco José de Caldas, Bogotá -Colombia. ⋆ ⋆ ⋆ Jannan David Rubio-Espinosa is undergraduate student of the Industrial Engineering Dept. of the Universidad Distrital Francisco José de Caldas, Bogotá -Colombia.


Hong-Quan Hua

T Hoang

Ngan Le

Bac Le

An Effective Initialization for ASM-based Methods

Keywords: Active Shape Models, Facial feature points, Initialization errors

Locating facial feature points is an important step for many facial image analysis tasks. Over the past few years, Active Shape Model (ASM) has become one of the most popular approaches to solve this problem. However, ASM-based methods are sensitive to initialization errors caused by poor face detection results. In this paper, an effective initialization for the ASM-based methods is proposed by our improved initial 8-landmarks ASM model together with Viola-Jones eye detection. In particular, we apply the 8-landmarks ASM model with the position of eyes from eye detector as reference points. After that, we choose several 76-landmark candidates from the training set that have the key feature points related to the result of the previous 8-landmark model. The best candidate has the lowest fitting errors with the test image and is used as initialization. To evaluate the performance of our work, we conduct the experiments on the MUCT and LFPW database. Compared to the latest ASM implementation, MASM, our proposed can improved the MASM by an average accuracy of 44% when dealing with poor initialization.

Introduction

Facial feature points play an important role for face analysis in many computer vision systems such as human-machine interface, face recognition and expression analysis. For this reason, the problem of facial feature points localization has attracted large amount of attention. As shown in the recent research by Chen, et al.'s [START_REF] Chen | Blessing of dimensionality: High-dimensional feature and its efficient compression for face verification[END_REF], leading performance in face verification can be acquired with simple features if facial feature points are detected accurately; therefore the problem of feature localization has become significant and attracted in the past years. Many approaches to extract facial feature points has been studied extensively including the active contour [START_REF] Kass | Snakes: Active contour models[END_REF], deformable template [START_REF] Yuille | Deformable templates for face recognition[END_REF], Active Shape Model [START_REF] Cootes | Active shape models-their training and application[END_REF], Active Appearance Models [START_REF] Cootes | Active appearance models[END_REF]. Most of above methods combine the local texture around each point and the shape of the whole face. ASM is known as one of the most popular methods that finds the landmarking points by Principle Component Analysis (PCA).

Active Shape Models (ASM) [START_REF] Cootes | Active shape models-their training and application[END_REF] proposed by Cootes, et al. has achieved great success in extracting facial feature points. In ASM, the feature points are modeled by the gray-scale texture whereas the geometrical information is defined by PCA. In order to implement a ASM system, a face detection is first used to find face regions and initialize a starting shape, then ASM is used to refine feature positions to the best locations based on the learned texture model. ASM is like deformable methods like snakes [START_REF] Kass | Snakes: Active contour models[END_REF], but ASM's feature points are constrained to be a valid shape by the learned statistical shape model. Though many methods have been proposed to improve the fitting accuracy of ASMs [START_REF] Cristinacce | Boosted regression active shape models[END_REF], [START_REF] Romdhani | A multi-view nonlinear active shape model using kernel pca[END_REF], [START_REF] Gu | A generative shape regularization model for robust face alignment[END_REF], facial feature points localization remains a very challenging problem. The challenge comes from the variations of facial appearance in illumination, pose, expression, aging and so on. Especially, ASM methods are very sensitive to the initialization perturbations (translation, rotation and scaling in size) caused by the poor face detection [START_REF] Seshadri | An analysis of the sensitivity of active shape models to initialization when applied to automatic facial landmarking[END_REF]. To overcome these defects, Cootes et al. [START_REF] Cootes | Multi-resolution search with active shape models[END_REF] propoesd a multi-resolution method, Wang et al. [START_REF] Wang | An improved active shape model for face alignment[END_REF] added edge information into the model. Milborrow and Nicolls [START_REF] Milborrow | Locating facial features with an extended active shape model[END_REF] achieved fairly accurate fitting results by combining 1D and 2D search model in the multi-resolution approach. Seshadri and Savvides [START_REF] Seshadri | Robust modified active shape model for automatic facial landmark annotation of frontal faces[END_REF] proposed a new robust method, called Modified Active Shape Model (MASM). MASM contains many theoretical changes to the traditional ASM and gains better results than above methods. However, MASM is still sensitive to initialization errors, especially while dealing with in-plane rotations and scaling effects. This is because the faces used to train the methods were upright and the initialization is too far from the optimal positions.

In our work, we focus on a new method to overcome the poor initialization. We select the initialization based on trained shapes instead of using mean shape as starting shape. Our selection process takes advantage of eye detector, support vector machine. Our selected feature points are almost in the optimal positions, which can prevent problems due to local minima. We also inherit the optimization metric of MASM [START_REF] Seshadri | Robust modified active shape model for automatic facial landmark annotation of frontal faces[END_REF] to determine the best location for each feature point. However, we weaken the important of edge information while fitting the facial boundary.

The remainder of this paper is organized as follows. In Section 2, the basic of Active Shape Model and improvements in Modified Active Shape Model are introduced. Our proposed method is described in Section 3. Experimental results are given in Section 4. A conclusion and future work are drawn in Section 5.

Related Works

This section gives an overview of Active Shape Models and the latest ASM-based methods called MASM. First, we explain the basic of training and testing stage in ASM method. Then, we show improvements of MASM to the basic ASM model since we use it as our main ASM-based method.

Active Shape Model

The following gives a brief review of the original ASM containing two main procedures: Training Stage and Testing Stage.

Training stage Shape Model and Local Appearance Model are built from a training set of manually annotated images. The Local Appearance Model is used to optimize the location for each landmark. On the other hand, the Shape Model constrains the landmark points to be in a valid shape.

The Shape Model : A landmark represents a distinguishable feature point in facial images. Locating landmarks is the same as locating feature points. Each face shape can be described by the coordinates of N landmarks points. The landmark points are (manually) annotated by a set of K training images. In our experiment, we use MUCT face database [START_REF] Milborrow | The muct landmarked face database[END_REF] which contains a total of 3755 faces with 76 manual landmarks at the resolution of 640 x 480 pixels. Figure 1 shows some example images of MUCT database. The coordinates of each landmarks for each image are stacked as a shape vector x.

x = [x 1 x 2 ...x N y 1 y 2 ...y N ] T

where x i and y i are the x and y coordinates of the ith landmark and N is the number of landmarks used.

The training images shape are all aligned by using Generalized Procrustes Analysis (GPA) [START_REF] Gower | Generalized procrustes analysis[END_REF] for minimizing the sum of squared distances between corresponding landmark points. PCA is then applied to reduce the dimensionality of the shape vectors, such that any shape x can be approximated as follows.

x ≈ x + Pb (2)
Here x is a mean shape vector, P is a matrix of eigenvectors and b = (b 1 b 2 ...b t ) T is a vector of shape parameters. When fitting the shape model to a face, the value of b i is constrained to lie with the range ±3 standard deviations to ensure that generated face shapes are valid.

The Local Appearance Model : The statistical models of appearance models (grayscale pixel intensities) of local regions around each landmark is used to optimize the location for each landmark in the search process. The 1D profiles are constructed by sampling the pixel intensities along a line, which is perpendicular to the landmark contour and centered at that landmark. The length of this line is 2k + 1 with k pixels sampled on either side of the landmark. The normalized derivative of this vector is computed to form the profile vector, g. The mean of such profiles (over the training images) at each landmark is called the mean profile vector, ḡ and the covariance matrix of such vectors is S g . To improve the accuracy and robustness, the process is repeated in a multi-resolution framework. The pyramid image contains L max levels, where image at level i is formed by smoothing the image at level i -1 and scaling by 50 percent. The fitting process is carried out at the coarse level first, then the results is refined in the finer level.

To compare the difference between a new sampled profile and the mean profile, we use the Mahalanobis distance measure as follows:

f (g i ) = (g i -ḡ) T S -1 g (g i -ḡ) (3)
By minimizing f (g i ), we can get the optimal position for the landmark. At each landmark position, a 2D profile is constructed by sampling a local patch around the landmark. The center of this local patch is the current landmark. The 2D local patch is then vectorized and normalized into vector g using a constant parameter c at each element of the profile, g i as shown in [START_REF] Cootes | Active appearance models[END_REF].

g i = g i g i + c (4)
where g i is the ith element of the vector g Testing Stage The search process of ASM can be implemented two iterative steps: (1) update the best location for each landmark point, (2) update the shape parameter b to ensure a valid shape. In the first step, we need to sample m pixels on both side of each landmark (m > k). We then compute the Mahalanobis distance as above for each sub-profile with the length of k. The sub-profile with the minimum distance is chosen and the the center of that sub-profile is the update position for the current landmark. Once all landmarks have been moved to the best position, we obtain the new shape vector x ini by using the shape model. The parameter vector b is calculated by minimizing the cost function in ( 5)

|x ini -T (x + Pb)| 2 (5)
In ( 5), b is the vector of shape parameter, T is the similarity transform that minimizes the Euclidean distance between x ini and x + Pb. The detail of the process to find b and T can be found in [START_REF] Cootes | Statistical models of appearance for computer vision[END_REF]. The adjusting process is repeated until no significant change in the landmark points in observed. Once convergence is reached at the lowest resolution level, the landmark coordinates is scaled and the entire process is repeated in lower level until convergence at the finest level. Figure 2 gives overview of ASM steps at testing stage.

Modified Active Shape Model

In this section, we describe the Modified Active Shape Model (MASM) proposed by Keshav Seshadri and Marios Savvides [START_REF] Seshadri | Robust modified active shape model for automatic facial landmark annotation of frontal faces[END_REF]. This approach gives significant improvements to the traditional ASM. MASM can handles faces with slight pose variations, in-plane rotations and varied expressions.

Firstly, MASM uses an optimal number of landmark points to represent face shapes. The landmarking scheme used in MASM contains 79 points which can perform accurately and model varied facial expressions. Secondly, edge detection is used to better fit points around the facial boundary. The main contribution of MASM is the development of new metric to determine the best location of the landmarks compared to the minimum Mahalanobis distance in traditional method.

MASM metric can be described as follows. MASM extracts 2D profiles and builds a subspace to model the pixel variation in the local patch around each landmark. They obtains the mean profile vector for each landmark and computes the eigenvectors of the covariance matrix S g corresponding to the dominant eigenvalues that models 97% of variation. Eigenvectors are stored along the columns of matrix P g . In the testing stage, a profile g around a candidate landmark is projected onto the subspace using equation ( 6) to obtain the shape parameter vector b g and is then reconstructed using this parameter vector to obtain g r using equation [START_REF] Cristinacce | Boosted regression active shape models[END_REF].

b g = P g T (g -ḡ) (6)
g r = ḡ + P g b g (7)
The metric for determining the best candidate point is the Mahalanobis distance between the original profile g and the reconstructed profile g r . The candidate with lowest reconstruction error is the best fit because it is the one that can be best modeled by the subspace. The reconstruction error can be computed as in equation [START_REF] Gower | Generalized procrustes analysis[END_REF].

error(g, g r ) = (g r -g) T S -1 g (g r -g) (8)
MASM also use edge information to improve the fitting accuracy for facial boundary. They assume that landmark points along the facial boundary usually lie along strong edges. The error function for these landmark is now defined in equation [START_REF] Gu | A generative shape regularization model for robust face alignment[END_REF].

error(g, g r ) = (c -I)(g r -g) T S -1 g (g r -g) (9)
where c is a constant, chosen to be 2 in MASM implementation and I is the Sobel edge intensity at the candidate point which can only take on the values 1 or 0. These modifications allow MASM to gain a significant improvement over the traditional ASM model.

Our proposed method

In this section, we introduce our proposed method for choosing initialization. In our method, each landmark is lie in an optimal location even when slightly expression and initialization errors are present. Moreover, we also introduce a modification to the cost function of MASM to improve the fitness of facial boundary. Figure 3 shows the comparison of initialization between our proposed method and MASM.

New method for choosing initialization

Most of previous ASM methods can not recover from poor initialization. The initialization is usually provided by the bounding box around the face which is the result of face detector. However, face detection algorithms such as Viola-Jones face detector [START_REF] Viola | Robust real-time face detection[END_REF] have no guarantees that the localization is perfect. There are many errors in translation and size. Especially, the bounding box from face detector is usually upright so the localization is poor when dealing with in-plane rotation of a face.

To overcome this problem, we propose a new method to reduce the sensitivity of ASM model to poor initialization. In our method, we use the location of eyes, nose and mouth to find 20 shape candidates from the training test that have the geometrical structure similar with the testing image. After that, we will use the local appearance information around each landmark to decide the best shape candidate. This shape candidate will be used as initialization instead of mean shape in previous ASM methods. Since we use the eyes, nose and mouth locations as well as the local appearance information, the resulting initial landmarks are almost in their optimal locations. The performance of this approach is shown in Section 4.

There are two different ASM model in our method, full-shape model and minimal-shape model. Full-shape model contains all the landmark points as in training set. On the other hand, minimal-shape model only contains 8 landmark points including 4 for eyes corners, 2 for nose and 2 for mouth corners. Figure 4 shows examples of both models. A flowchart of the overall process is depicted in Fig. 4: Examples of full-shape model on the left and mini-shape model on the right Figure 5. First, we use Viola-Jones face detector to extract the face regions from the test image. Then, we detect the eyes using Viola-Jones algorithm. Since the detected regions contain many false positives, we apply a simple linear Support Vector Machine (SVM) to select true positive eye regions. After that, we apply the minimal-shape model with eyes positions as reference points. The result is the optimal positions for 8 landmark points and the corresponding shape parameter vector b min . The key idea is that we will use this b min to select 10 nearest shape Fig. 5: Sequence of operations to select the suitable initialization for MASM model parameter vector of known minimal shape in training set. This means that we now have 10 candidates of full shape for initialization purpose.

Because the eye detector doesn't guarantee perfect results, we repeat the above process to have 10 more candidates of full shape. However, we use fixed positions in the bounding box as reference points for minimal-shape model instead of eyes positions. The ratio for fixed positions is described in Figure 6. After obtaining 20 candidates for initialization, we will select the candidate that has lowest error value. The error value of each candidate is the summation of fitting error of the profile around each landmark point. The selected candidate will be use as initialization points for full-shape model.

Use of edge information

In MASM, Keshav Seshadri and Marios Savvides assume that facial landmark usually lie along strong edges. However, we suggest that we should weaken that assumption because Sobel edge detector don't guarantee that strong edges will lie on facial boundary. In that case, the candidate points may be chosen inaccurately and then the performance of fitting process is reduced.

For the landmarks in the facial boundary, we will use the cost function as in equation [START_REF] Kass | Snakes: Active contour models[END_REF] to measure Mahalanobis distance between the original profile g and the reconstructed profile g r .

error(g, g r ) = (1 -c * I)(g r -g) T S -1 g (g r -g) (10)
where c is a constant, chosen to be 0.2 in our implementation, I is the Sobel edge intensity at the candidate point and S g is the covariance matrix of the current landmark's subspace.

Experiments and Results

In this section, we will evaluate the performance of our proposed method and show the experimental results. The comparison experiments between our method and the Modified Active Shape Model (MASM) are conducted on MUCT [START_REF] Milborrow | The muct landmarked face database[END_REF] and LFPW [START_REF] Belhumeur | Localizing parts of faces using a consensus of exemplars[END_REF] database. All model are trained with 300 face images from MUCT database. The landmark scheme used in MUCT database contains 76 landmark points as in Figure 4. Some original sample images of MUCT database are shown in Figure 1. We conducted two experiments in testing stage. In the first test, we evaluate with two test set from MUCT and LFPW with the assumption that initialization errors were not present. The first set consisted of 450 images from the MUCT database. The second set consisted of 300 images from LFPW database. In the second experiment, we added various of initialization errors to the face detector results.

To evaluate the performance of the method, we use the distance between the predicted landmark positions and the manually annotated points.

E = 1 n n i=1 { 1 N N j=1 dist(s ij , m ij )} ( 11
)
where n is the number of test images, N is the number of landmark, s ij and m ij are the coordinate vectors of the automatic landmark and the manually labeled landmark respectively. The Euclidean distance dist(s ij , m ij ) between s ij and m ij is normalized so that the distance between the nearest corners of both eyes is 30 pixels. For these experiment, we used Viola-Jones face detector to get the face bounding box. Then, we refined manually the incorrectly detected results in translation and scaling to prevent initialization errors. The results of first test are illustrated in Table 1. We can find that our method can enhance the accuracy of the MASM by 23.66 % in MUCT test and 14.59% in LFPW test. Figure 7 shows plots the errors for each landmark. The landmark scheme of MUCT database is 76 points. However, the ground truth landmark points in LFPW is 68 points, as labeled by IBUG [START_REF] Sagonas | A semi-automatic methodology for facial landmark annotation[END_REF]. Therefore, we choose the corresponding landmark points in both scheme to compute fitting errors. The resulting scheme contains 61 points. As can be seen, our proposed method has better result than the MASM method, especially while fitting the landmarks along the facial boundary [START_REF] Belhumeur | Localizing parts of faces using a consensus of exemplars[END_REF][START_REF] Chen | Blessing of dimensionality: High-dimensional feature and its efficient compression for face verification[END_REF][START_REF] Cootes | Multi-resolution search with active shape models[END_REF][START_REF] Cootes | Active appearance models[END_REF][START_REF] Cootes | Active shape models-their training and application[END_REF][START_REF] Cootes | Statistical models of appearance for computer vision[END_REF][START_REF] Cristinacce | Boosted regression active shape models[END_REF][START_REF] Gower | Generalized procrustes analysis[END_REF][START_REF] Gu | A generative shape regularization model for robust face alignment[END_REF][START_REF] Kass | Snakes: Active contour models[END_REF][START_REF] Milborrow | The muct landmarked face database[END_REF][START_REF] Milborrow | Locating facial features with an extended active shape model[END_REF][START_REF] Romdhani | A multi-view nonlinear active shape model using kernel pca[END_REF][START_REF] Sagonas | A semi-automatic methodology for facial landmark annotation[END_REF][START_REF] Seshadri | Robust modified active shape model for automatic facial landmark annotation of frontal faces[END_REF]. This is the effect of the modification in using edge information. On the remaining landmark points, our method just gives slightly improvements because we focus on the problem of initialization errors while these experiments assume that there is no initialization errors. Besides, other processing steps are almost the same as MASM. Hence, we don't have much improvement on others part. In the LFPW test, the errors are higher than the first test because LFPW images are captured in real-world scenario which contain many pose changes and expressions. However, our method can gain better results than MASM.

For our second experiments, we estimate the performance of our method when dealing with initialization errors. In this experiments, we test the effect of translation and scaling to the accuracy of our proposed method and MASM. Figure 8 shows the results of our test. In the translation test, we move the face bounding box with a same value along x and y coordinates. The translation value is: -50, -25, 0, 25, 50. In the scaling test, we scale the face bounding box with ratios of [ 1.5, 1.2, 1, 0.8, 0.6 ]. In Scaling test, we achieve the fitting errors of 2.32048 compared to 3.9663 of the MASM. In translation test, the error of our method is 1.8742 while MASM's error is 3.4945. In both cases, our proposed method can gains better performance, 46.37% and 41.49% in translation and scaling test respectively , which means that our method are more robust to the initialization errors.

Conclusion and Future work

In this paper, we presented a novel Active Shape Model for facial feature points detection in frontal-view face images. Our method focus on dealing with initialization errors due to the failure of face detector. In order to recover from poor initialization, we use an eye detector with Support Vector Machine to effectively extraction the location of the eye. Then, we apply a minimal ASM model with 8 key landmark points to obtain an optimal initialization. This optimal initialization is used as the starting point for the full ASM model with 76 landmark points instead of using mean shape. Moreover, we consider the usage of edge information carefully to prevent misalignment on facial boundary. Experimental results show that our proposed method can achieve better fitting accuracy even with poor initialization. In the future work, there are still many rooms to improve the performance of our proposed method. First, our method uses the eye location to obtain the initialization. Hence, our method can't deal with occlusion effectively, especially when the target face has sun glasses. Such problems can be solve by combing the location of nose and mouth in addition to eye location. Secondly, the feature descriptor is very simple. We can improve it with HOG or Gabor features which are robust to illumination.

Fig. 1 :

 1 Fig. 1: Example images from MUCT database. The green dots indicate the manually labeled landmark points

Fig. 2 :

 2 Fig. 2: Overview of ASM steps at testing stage. The blue box is the result of face detector. Red dots indicate the starting shape for each step and Green dots indicate the landmark result

Fig. 3 :

 3 Fig. 3: Initialization of MASM and our proposed method without initialization errors a) and with initialization errors b). The Red dots indicate the initialization of MASM. The Green dots indicate the initialization of our proposed method

Fig. 6 :

 6 Fig. 6: The position of fixed minimal-shape points in the face region

Fig. 7 :

 7 Fig. 7: Plots of fitting errors vs landmark number for (a) MUCT test (b) LFPW test.

Fig. 8 :

 8 Fig. 8: Plots of fitting errors with the effects of scaling a) and translation b)

Table 1 :

 1 Results for first test without initialization errors

		Fitting Errors Our method	MASM	Improvement
	MUCT	1.5777	2.0666	23.66%
	LFPW	7.6756	8.9873	14.59%
	(a)			(b)

Acknowledgments

This research is supported by research funding from Advanced Program in Computer Science, University of Science, Vietnam National University -Ho Chi Minh City.