R Latha
email: lathapsg@yahoo.co.in

R Nadarajan
email: nadarajan_psg@yahoo.co.in

User Relevance for Item-based Collaborative Filtering

Keywords: Collaborative filtering, Recommendation System, Information Retrieval, User Relevance

A Collaborative filtering (CF), one of the successful recommendation approaches, makes use of history of user preferences in order to make predictions. Common drawback found in most of the approaches available in the literature is that all users are treated equally. i.e., all users have same importance. But in the real scenario, there are users who rate items, which have similar rating pattern. On the other hand, some users provide diversified ratings. We assign relevance scores to users based on their rating pattern in order to improve the quality of predictions. To do so, we incorporate probability based user relevance scores into the similarity calculations. The improvement of predictions of benchmark item based CF approach with the inclusion of user relevance score is demonstrated in the paper.

Introduction

With the tremendous growth of Web, volume of data available in net based systems become large and thus results in large repositories. So it becomes too difficult for individuals to handle the data effectively and efficiently. This scenario is commonly referred to as 'information overload' problem. Recommender system addresses the problem [START_REF] Melville | Recommender Systems, Encyclopedia of Machine Learning[END_REF]. Many online business systems such as Amazon.com and Netflix are using recommender systems to provide personalized suggestions.

Two common approaches are available for making recommendations. They are Collaborative Filtering (CF) [START_REF] Breese | Empirical Analysis of Predictive Algorithms for Collaborative Filtering[END_REF] and Content Based filtering (CB) [START_REF] Pazzani | Content-based recommendation systems: The Adaptive Web[END_REF]. Pazzani, et. al [START_REF] Pazzani | Content-based recommendation systems: The Adaptive Web[END_REF] defines that content-based algorithms base their recommendations on the contents of items and profiles of users. The profiles allow programs to associate users with matching products. Content Based Filtering uses the assumption that items with similar objective features will be rated similarly and users with similar profile/taste will prefer items in a similar manner. CB approaches are specific to a domain and the scope of the approach is limited to the domain for which they are proposed.

Collaborative Filtering (CF) is the most popular and successful approach for recommendation systems. CF relies only on the past user behavior (ratings, preferences, purchase history, time spent etc) [START_REF] Melville | Recommender Systems, Encyclopedia of Machine Learning[END_REF]. Breese et. al [START_REF] Breese | Empirical Analysis of Predictive Algorithms for Collaborative Filtering[END_REF] classifies Collaborative Filtering techniques into two categories: Memory-based and Model based techniques. Memory based algorithms predict the rating of users using the previously rated items by the users and other users who have similar tastes. They operate over the entire user database to make predictions. The most common memory-based models are based on the notion of nearest neighbors, using a variety of distance measures [START_REF] Pazzani | Content-based recommendation systems: The Adaptive Web[END_REF]. Model based systems are based on a compact model inferred from the data. We have considered model based CF approach for our work.

In Collaborative Filtering techniques, all users are treated uniformly. But in the real scenario, there are some users who are consistent in providing preferences, where as some users provide ad hoc preferences which do not form any pattern. Ultimately the former kind of users must be assigned higher score, and so their preferences must be given higher importance. We assign scores to users based on their rating pattern in order to improve the quality of predictions. We assign probability based relevance score to users.

The contribution of the proposed work can be divided into two phases, namely offline phase and online phase. In offline phase user relevance scores are calculated and a model is built on relevance scores. In on line phase rating predictions are made based on the model created.

The rest of this paper is organized as follows.

Section 2 presents the overview of item based CF and user relevance problem, Section 3 discusses about the proposed approach to address the problem. Section 4 presents user relevance model building and Section 5 discusses about experimental evaluations of our approach. MovieLens [START_REF]MovieLens data[END_REF] database is used for proving the results. Finally, Section 6 provides some concluding remarks and an outline of the future research.

Overview of Item Based CF and User Relevance Problem

In this section we discuss about Item based CF techniques and user relevance problem.

Item Based CF

To address the scalability concerns of user-based recommendation techniques, item-based recommendation techniques (also known as model-based) have been developed [START_REF] Karypis | Evaluation of Item-Based Top-N Recommendation Algorithms[END_REF]. These techniques compute similarity between items, and then use these similarity values to compute top-N recommendations or make predictions. The reason behind these techniques is that a customer will more likely purchase items that are similar or related to the items that he/she has already purchased. These approaches are faster than user based approaches, since the similarity computations can be done offline which results in faster recommendation engines.

Different kinds of similarity measures are available in the literature for computing relationship between items. In [START_REF] Sarwar | Item-based collaborative filtering recommendation algorithms[END_REF], the author has discussed about cosine similarity, adjusted cosine similarity and Pearson correlation coefficient.

In Cosine based approach [START_REF] Adomavicius | Toward the Next Generation of Recommender Systems: A Survey of the State-of-the-Art and Possible Extensions[END_REF] the similarity between items i and j is calculated as cosine of the angle between them. Only the items which are rated by both the users will be considered for computing the angles. The similarity is defined as

∑ √ ∑ √ ∑
where represents set of users who have rated for the item i and represents the ratings given by the user u for the item i. According to B. Sarwar et.al, [START_REF] Sarwar | Item-based collaborative filtering recommendation algorithms[END_REF] Cosine similarity does not account for the difference in user ratings. Adjusted Cosine Similarity is suggested to overcome the problem, in which the average rating of the user is subtracted from his actual ratings. The function used to calculate the similarity between item i and item j is given below:

∑ ( ̅ ) ̅ √ ∑ ̅ √ ∑ ̅
where represents the users who have rated for the item i and represents the ratings given by user u for item i. ̅ is the average rating of user u. There are many ways to compute predictions [START_REF] Adomavicius | Toward the Next Generation of Recommender Systems: A Survey of the State-of-the-Art and Possible Extensions[END_REF]. Weighted sum is a commonly used approach for predicting unknown ratings which is calculated as given in (3)

∑ ∑ | |
Here S(i) is the set of items that are similar to i, the item to be predicted. is the prediction for item i for the user u. Although CF techniques take different kinds of input, we focus on the input in the form of m x n user item matrix.

Definition 1. User-item Matrix R

If there are m users who have given ratings for n items, then the ratings data can be represented as an m x n matrix with rows representing users and columns representing items. The matrix is called user-item matrix R. Each element is an ordinal value which ranges from to . Unrated values are considered to be zero. A sample rating matrix R is shown in Table 1.

User Relevance in Item Based CF

In Collaborative Filtering techniques we discussed above, the importance of users is not taken into account. There are some users who prefer items with similar rating pattern. i.e, if a user rates an item and most of its similar items, then the user is consistent in rating. By assigning higher weight to them brings accuracy in recommendations. On the other hand there are some users who rate items which do not have any similar rating pattern. They rate dissimilar items. Such users are divergent in thinking and assigning higher weight to them improves diversity in recommendations. But current recommender systems treat all users uniformly. Min Gao et. al [START_REF] Gao | Userrank for item-based Collaborative filtering recommendation[END_REF] proposed a technique to assign rank to users based on popularity and use the rank values in similarity calculations. The technique improves accuracy of predictions.

Our approach proposed in this paper is similar to the work described in [START_REF] Gao | Userrank for item-based Collaborative filtering recommendation[END_REF], but instead of user rank, user relevance score is calculated. Moreover the proposed technique improves accuracy as well as diversity of recommendations.

Despite of much work available in the literature to improve quality of predictions, up to our knowledge no other work is available in the literature to assign weights to users and treat them differently.

Table 1. Original rating Matrix R I1 I2 I3 I4 I5 I6 U1 1 2 0 5 0 0 U2 0 3 0 0 5 0 U3 0 0 5 0 0 2 U4 0 1 0 0 0 0

Probability based Relevance Score for Users

This section describes the procedure for building a model for assigning relevance score for the users. To calculate the relevance score for users we apply a probability based technique. The model building takes two phases namely item-item similarity calculation phase and relevance score computation phase.

Item-item Similarity Calculation Phase

In order to define pair-wise similarity between item vectors many similarity measures namely Pearson correlation coefficient, cosine similarity, adjusted cosine similarity etc are available in the literature. The similarity measure we have considered is based on number of users who agree on item vectors i and j. The Correlation Rating Matrix (CRM) stores the similarities between each pair of item vectors. between item vectors i and j is defined as given in (4)

{ | |
where ui is the set of users who have rated for item i. | | represents cardinality of the set ( )

Definition 2. Correlation Rating Matrix (CRM)

For the given m × n user item rating matrix R, item-item Correlation Rating Matrix can be represented as an n × n matrix. The matrix rows and columns represent items and each represents 1 if at least one user rates both the items. = 0, implies that item i is not related to item j.

CRM is a symmetric matrix. But if an item i is rated by few users and item j is rated by many users, then their correlation should differ. So we normalize CRM in to Normalized Correlation Similarity Matrix NCRM.

Definition 3. Normalized Correlation Rating Matrix (NCRM)

NCRM is a matrix that records the relationship between items as ratio of how the user is correlating in the current item and the remaining items. NCRM can be calculated based on the formula given in (5)

∑

The rows and columns of NCRM represent items and it is an asymmetrical matrix. For each item, the NCRM value is between 0 and 1. Since the relevance score computation procedure we employ requires binary similarity matrix for relevance score computations. So we convert NCRM into S, a binary similarity matrix of various items of the rating matrix. CRM and NCRM of original rating matrix R are shown in Fig. 1.

Fig. 1. CRM and NCRM of original matrix R given in Table1 Definition 4. Item-item Similarity Matrix (S)

For the given m × n user item rating matrix R, item-item similarity matrix is a binary similarity matrix(S) which can be calculated by the formula { S is a binary similarity matrix. Similarity matrix S of R is shown in Table 2.

Table 2. Binary Similarity Matrix, S I1 I2 I3 I4 I5 I6 I1 0 0 0 0 0 0 I2 0 0 0 0 0 0 I3 0 0 1 0 0 1 I4 0 0 0 0 0 0 I5 0 1 0 0 1 0 I6 0 0 1 0 0 1 3.2
Relevance Score Computation Phase.

From the binary similarity matrix S, for each item, the items with similarity value 1 are called relevance items (R) and with similarity value 0 are called non-relevant items (NR). Next we investigate how user u has rated an item i and all its similar and dissimilar items. If he has rated most of the similar items of i, then he can be assigned a good score for his consistency in rating. At the same time high consistence reveals monotonic rating behavior of users which needs to be balanced.

On the other hand if the user has rated an item and most of the dissimilar items of it, then he can be assigned a good score for his divergent thinking and rating. So we consider both scores in order to have a tradeoff between them. We combine both of them to calculate relevance score of the user. The relevance score of the user is estimated using [START_REF] Ekstrand | Collaborative filtering Recommendation systems[END_REF] and [START_REF] Pazzani | Content-based recommendation systems: The Adaptive Web[END_REF]. This is based on the basic probability model proposed by Robertson and Spark [15] in Information Retrieval domain. This feature weighing scheme is used by search engines to rank matching documents according to their relevance to a given search query.

where and are the consistency and diversity scores of user u for item i. and are log odd ratios of user u in relevant and non-relevant items of item i respectively. N is the number of items in the collection. R is the number of items similar to item i. NR is the total number of dissimilar items of i.

is the num-ber of items in the collection for which the user u has rated and is the number of similar items of i for which user u has rated, n is the number of dissimilar items of i for which user u has rated. Weight assigned to each user is calculated as the total scores assigned to all items and is calculated as given in ( 9) and [START_REF] Su | A Survey of Collaborative Filtering Techniques[END_REF].

∑ ∑

where I is the set of items and and are called consistency and diversity scores of user u. Consistency score signifies user's consistency in rating items where as diversity signifies user's diversified thinking in analyzing items. Since both the scores are important as one improves accuracy of predictions where as other improves diversity in predictions, we consider average of the two scores as a common score for users, which we call as relevance score. Relevance score of users is calculated as given in [START_REF] Gao | Userrank for item-based Collaborative filtering recommendation[END_REF].

Building User Relevance Model for CF

In the previous section we discussed about computing user relevance score. We incorporate user relevance score into the most widely used similarity measure namely Cosine Similarity. A new similarity measure which we call User Relevance based Cosine Similarity, is formulated as given in [START_REF] Lemire | Slope One Predictors for Online Rating-Based Collaborative Filtering[END_REF].

∑ ( ) ( ) √∑ √∑

RWu is the relevance score of the user who rated both the items i and j. After calculating the similarity values, a model is built as described in [START_REF] Sarwar | Item-based collaborative filtering recommendation algorithms[END_REF].

Experimental Evaluation

This section discusses about the data sets used, accuracy of predictions and computational complexities. Experiments are conducted using Item-based collaborative filtering recommendation algorithms [START_REF] Sarwar | Item-based collaborative filtering recommendation algorithms[END_REF] and Userrank for item-based Collaborative filtering recommendation [START_REF] Gao | Userrank for item-based Collaborative filtering recommendation[END_REF], in order to prove the efficiency of the proposed technique.

Data Set Used

The experiments are run on two datasets namely Movielens [START_REF]MovieLens data[END_REF] and a subset of rating selected from it. The subset is populated by taking top 100 users based on total number of rating and top 100 items rated by those users. We call it as MovieTop100. The datasets have ratings given for movies in the range 1 to 5.

The details of the data sets are given in Table 3. We partition the data sets into training set with 80% of the ratings and test data set with 20% of the ratings. Five cross validation is done for all experiments. Here we outline the experiments we have taken to study the significance of the proposed technique in improving the quality of predictions in terms of accuracy and diversity.

Accuracy Measures

In order to evaluate the performance of the proposed technique in terms of accuracy of predictions, we follow the approach used in [START_REF] Yildirim | A Random Walk Method for Alleviating the Sparsity Problem in Collaborative Filtering, RecSys[END_REF]. The metrics used are HR (Hit Ratio), MAE (Mean Absolute Error) and RMSE (Root Mean Square Error). HR is the ratio of the number of hits to the size of test data set used. The predicted rating is a called as a hit if its rounded value is equal to the actual rating given in the test data set. MAE is average absolute deviation between predicted and actual ratings. MAE penalizes each wrong prediction by its distance to the actual rating, whereas RMSE emphasizes larger deviations. These measures are as formulated below.

∑ | | √ ∑

Where p i is the predicted rating and a i is the actual rating for the item i. n is the number of test cases.

Diversity Measure (Cosine diversity)

In order to show that the proposed approach improves diversity of predictions we employ average inter list dissimilarity of predicted values. The cosine coefficient is a commonly used similarity measure. Based on cosine similarity a new diversity measure called cosine dissimilarity is proposed in [START_REF] Hurley | Novelty and diversity in top-n recommendation -Analysis and evaluation[END_REF]. Cosine dissimilarity is calculated for each pair of predicted ratings and its average is taken as diversity value.

We compared our technique with ItemCF [START_REF] Sarwar | Item-based collaborative filtering recommendation algorithms[END_REF] and ItemCF(U-Rank) [START_REF] Gao | Userrank for item-based Collaborative filtering recommendation[END_REF]. Table 4 shows the results of computations for various measures. The quality of predictions increases as HR and Diversity increase whereas quality decreases as MAE and RMSE increase. The values of various measures from Table 4 show that the proposed approach performs better than the benchmark techniques ItemCF and ItemCF(U_Rank) on both the datasets. High computational complexity is often needed to enhance the predictions. In the model based view, the computational complexity can be split into complexity in offline phase and complexity in online phase. Offline phase includes similarity computations, relevance score calculations and model building. The complexity of offline computations is O(mn2+n2+kmn) for m users, n items and top k neighbours. Online phase includes prediction computations of n test items with complexity O(kn). For most recommender systems, the online complexity is more important than offline complexity.

Conclusion

The goal of this work is to improve accuracy and diversity of predictions by including relevance score of users. The proposed approach sets user relevance score as a combination of user's consistency and diversity scores and a new similarity measure is proposed based on user relevance score. A bench mark algorithm with the inclusion of user relevance score is examined for proving the efficiency of results. The experimental results show that the inclusion of relevance weight helps in improving predictions. In future we plan to include the content based attribute values in order to calculate the relevance score of the users.

Fig. 2 and

 2 Fig.2and Fig.3show the comparison of all the three algorithms using different measures on Movielens dataset.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Comparison of MAE and RMSE measures for Movielens dataset

Fig. 4 .

 4 Fig. 4. Comparison of MAE and RMSE measures for MovieTop100 dataset

Fig. 5 .

 5 Fig. 5. Comparison of HR and Diversity measures of MovieTop100 dataset

Table 3 .

 3 Data sets used

			Movielens	MovieTop100
		# users	943	100
		# items	1682	100
		# ratings	100,000	2474
	5.2	Experimental Metrics and Evaluation Methodologies.

Table 4 .

 4 Comparison of various measures by ItemCF, ItemCF(U_Rank) and ItemCF(U_Rel)

			Movielens			MovieTop100
		ItemCF	ItemCF	ItemCF	ItemCF	ItemCF	ItemCF
			(U_Rank)	(U_Rel)		(U-Rank)	(U_Rel)
	MAE	0.9318	0.9306	0.9155	0.9208	0.9055	0.9020
	RMSE	1.2041	1.1976	1.1778	1.1710	1.1575	1.1519
	HR	0.3386	0.3388	0.3441	0.3233	0.3278	0.3288
	DIV	0.0243	0.0244	0.0250	0.0152	0.0175	0.0177


Krzysztof Gdawiec
email: kgdawiec@ux2.math.us.edu.pl

Aesthetic Patterns from the Perturbed Orbits of Discrete Dynamical Systems

Keywords: dynamical system, orbit, perturbation mapping, aesthetic pattern

The aim of this paper is to present some modifications of the orbits generation algorithm of discrete dynamical systems. The first modification is based on introduction of a perturbation mapping in the standard Picard iteration used in the orbit generation algorithm. The perturbation mapping is used to alter the orbit during the iteration process. The second modification combines the standard Picard iteration with the iteration which uses the perturbation mapping. The obtained patterns have unrepeatable structure and aesthetic value. They can be used for instance as textile patterns, ceramics patters or can be used in jewellery design.

Introduction

One of the most elusive goals in computer aided design is artistic design and pattern generation. Pattern generation involves diverse aspects: analysis, creativity, development [START_REF] Wannarumon | Intelligent Computer System for Jewelry Design Support[END_REF]. A designer have to deal with all of these aspects in order to obtain an interesting pattern which later could be used in jewellery design, carpet design, as a texture etc. Usually the most work during the design stage is carried out by a designer manually. Especially, in the case in which a graphic design should contain some unique unrepeatable artistic features. Therefore, it is highly useful to develop an automatic method for aesthetic pattern generation.

Aesthetics in the world of art and photography is connected with the principles of the nature and the perception of beauty [START_REF] Datta | Studying Aesthetics in Photographic Images Using a Computational Approach[END_REF]. Judging the beauty and other aesthetic qualities of patterns, paintings, photographs is a highly subjective task, so there is no standard method of measuring aesthetic values and it is a challenge to create such method in the emerging discipline of computational aesthetic [START_REF] Wu | The Good, the Bad, and the Ugly: Predicting Aesthetic Image Labels[END_REF]. However, in most of the works about pattern generation, in this paper also, the aesthetic is judged by the subjective feeling of the authors.

The literature is full of very diverse methods of pattern generation. For instance in [START_REF] Pang | Interactive Evolutionary 3D Fractal Modeling[END_REF] [START_REF] Wannarumon | A New Aesthetic Evolutionary Approach for Jewelry Design[END_REF] the methods based on Iterated Function Systems and Genetic Algorithms for jewellery design were proposed. An interesting method based on root-finding polynomials, called polynomiography, was presented in [START_REF] Kalantari | Polynomial Root-Finding and Polynomiography[END_REF]. In [START_REF] Sen | A Product-Delay Algorithm for Graphic Design[END_REF] a product-delay algorithm was proposed which is based on multiplication of sine and square waves and in [START_REF] Lu | Mathematical Marbling[END_REF] authors presented a method which creates stone-like decorations using marbling.

Also very diverse methods based on discrete dynamical systems were proposed. In [START_REF] Naud | Automatic Generation of Aesthetic Images for Computer-assisted Virtual Fashion Design[END_REF] authors made use of Gumowski-Mira transform for the generation of patterns used later as textile patterns. In their method they took the Gumowski-Mira transform and changed only the parameters of it. Other approach was presented in [START_REF] Gdawiec | Automatic Generation of Aesthetic Patterns with the Use of Dynamical Systems[END_REF]. The authors have used different kinds of dynamical systems and using the switching process between the dynamical systems in conjunction with the Krasnoselskij iteration they generated aesthetic patterns. In this paper we propose another approach in obtaining aesthetic patterns from discrete dynamical systems. We modify the standard Picard iteration using the perturbation mapping and also combine both methods.

The paper is organized as follows. In Sect. 2 we introduce the basic information about discrete dynamical systems and present examples of dynamical systems which produce nicely looking orbits. Next, in Sect. 3 we introduce the notion of a perturbation mapping and how it can be used in the pattern generation. Some examples of aesthetic patterns obtained with the proposed methods are presented in Sect. 4. Finally, in Sect. 5 we give some concluding remarks.

Discrete Dynamical Systems

Let M be a subset in the q-dimensional Euclidean space IR q . Following [START_REF] Osipenko | Dynamical Systems, Graphs, and Algorithms[END_REF] discrete dynamical system is a continuous mapping Φ :

M × IN → M such that Φ(x, 0) = x, (1)
Φ(Φ(x, t), s) = Φ(x, t + s) (2)
for all t, s ∈ IN and x ∈ M . The variable t is thought of as the time and generally discrete dynamical systems result from iterative processes or difference equations [START_REF] Osipenko | Dynamical Systems, Graphs, and Algorithms[END_REF]. Assume that f : M → M is a continuous mapping. Then f generates a discrete dynamical system of the form [START_REF] Osipenko | Dynamical Systems, Graphs, and Algorithms[END_REF]:

Φ(x, n) = f n (x) = f • . . . • f n times (x). ( 3
) If n = 0 then f 0 (x) = x.
In the rest of the paper we will be interested in discrete dynamical systems given in a following form:

x n = f n (x) = f (f n-1 (x)) = f (x n-1 ), (4)
where n > 0. Formula ( 4) is also called the Picard iteration.

For our further considerations we also need the notion of an orbit. Orbit (or trajectory) of a point x 0 is a sequence {x n } ∞ n=0 , where x n is given by [START_REF] Kalantari | Polynomial Root-Finding and Polynomiography[END_REF]. Many examples of dynamical systems are known [START_REF] Morozov | Invariant Sets for Windows[END_REF], but we are mainly interested in those which produce geometric patterns that can be recognized as aesthetic ones. Now, we present the examples of such dynamical systems in IR 2 :

-Hopalong transformation [START_REF] Martin | Graphic Potential of Recursive Functions[END_REF] x n = y n-1 -sgn(x) |bx n -c|,

y n = a -x n-1 , (5)
where a, b, c ∈ IR and sgn : IR → IR is defined as follows:

sgn(x) =      -1 if x < 0, 0 if x = 0, 1 if x > 0, (6)
-Zaslavsky transformation [START_REF] Morozov | Invariant Sets for Windows[END_REF] x n = (x n-1 + K sin y n-1 ) cos α + y n-1 sin α,

y n = -(x n-1 + K sin y n-1 ) sin α + y n-1 cos α, (7)
where K ∈ IR, α = 2π q , q ∈ IN, q ≥ 3, -Chip transformation created by Peters for the HOP program [START_REF] Peters | HOP -Fractals in Motion[END_REF] x

n = y n-1 -sgn(x n-1 ) cos(ln |bx n-1 -c|) 2 • arctan(ln |cx n-1 -b|) 2 , y n = a -x n-1 , (8)
where a, b, c ∈ IR, -Quadrup Two transformation also created by Peters for the HOP program [START_REF] Peters | HOP -Fractals in Motion[END_REF] x n = y n-1 -sgn(x n-1 ) sin(ln

|bx n-1 -c|) arctan(cx n-1 -b) 2 , y n = a -x n-1 , (9)
where a, b, c ∈ IR, -Three Ply transformation is another transformation created by Peters [START_REF] Peters | HOP -Fractals in Motion[END_REF] x

n = y n-1 -sgn(x n-1 )| sin x n-1 cos b + c -x n-1 sin(a + b + c))|, y n = a -x n-1 , (10)
where a, b, c ∈ IR.

The other examples of dynamical systems, producing interesting orbits, are: Gumowski-Mira transformation (CERN, 1980) [START_REF] Gumowski | Recurrences and Discrete Dynamic Systems[END_REF], Martin [START_REF] Martin | Graphic Potential of Recursive Functions[END_REF], Cockatoo [START_REF] Martin | Graphic Potential of Recursive Functions[END_REF] etc. In Fig. 1 the examples of orbits for the first 100000 points for transformations (5), ( 7), ( 8), ( 9), [START_REF] Pang | Interactive Evolutionary 3D Fractal Modeling[END_REF]

Orbit Perturbation

As we mention at the beginning of the paper there exist methods which use discrete dynamical systems presented in Sect. 2. They use changing of the transformation parameters or switching process in conjunction with the Krasnoselskij iteration. In our approach we are interested in the orbits of discrete dynamical systems.

Lets see an example of orbits for different starting points of a fixed dynamical system. Figure 2 presents orbits starting from (from left): (1.2, -2.05), (-3.55, -1.2) and (3.601, -1.05) for the Quadrup Two transformation with a = 34, b = 1, c = 5 and n = 100000. To colour the orbits we used the iteration method [START_REF] Gdawiec | Automatic Generation of Aesthetic Patterns with the Use of Dynamical Systems[END_REF] in which the colour of a point is given according to the number of iteration at which this point arises and a fixed colourmap. In this example we see that for different starting points we obtain various shapes of the orbits. This gives an idea that we can alter the orbits during the iteration process to change the shape of the final orbit.

We modify the Picard iteration given by (4) in a following way:

x n = (f • p) n (x) = (f • p)((f • p) n-1 (x)) = (f • p)(x n-1 ), (11)
where p : M → M is a mapping. The mapping p is called a perturbation mapping and its aim is to alter (perturb) the orbit during the iteration process. We do not make any assumptions about the perturbation mapping because we can alter A very simple example of perturbation mapping is:

p v (x) = x + v, (12)
where v ∈ M . This mapping simply translates the given point by v. In this way we change the orbit in every iteration in the same direction. We also can change the orbit by v only in some iterations, e.g. if the number of iteration is divisible by some fixed value.

Of course there is nothing in the way to use different values of v in consecutive iterations. In this case we can define some equation for v or we can define a grid with vectors and the appropriate vector is then computed using the bilinear interpolation. We also can give some randomness in the perturbation mapping taking in each iteration a random value of v.

When we have a perturbation mapping p we can define a new iteration process which combines the standard Picard iteration with the iteration given by [START_REF] Peters | HOP -Fractals in Motion[END_REF]. For instance this can be done in a following way:

x n = αx n + (1 -α)x n , (13)
where α ∈ IR and

x n = f (x n-1 ), (14)
x n = f (p(x n-1 )). ( 15
)

Examples

In this section we show some examples obtained using the methods presented in Sect. [START_REF] Gumowski | Recurrences and Discrete Dynamic Systems[END_REF]. In all the examples we use the same colouring method as in Fig. 2 but with other colourmaps. We start with an example of Zaslavsky transformation with K = 3, q = 5, starting point (1, 1) and number of iterations equal 100000. In the iteration process we used p v given by ( 12) with v = (1, 0.5) which is used after every: 500, 2500, 10000 iterations. The obtained results are presented in Fig. 3. The top part of the figure presents the original orbit without any perturbation, and the bottom part presents the patterns obtained with the perturbed orbits after every (from left): 500, 2500, 10000 iterations. Using the random vectors in the perturbation mapping we are also able to obtain very interesting and diverse patterns which shows next example. For this example we used the Hopalong transformation with a = -55, b = -1, c = -42, starting point (0, 0) and number of iterations equal 100000. The obtained patterns are presented in Fig. 5, where the original orbit is in the top part, and in the bottom part we have patterns obtained with the random vectors. The co-ordinates of the random vector were randomly chosen from [-0.05, 0.05].

In the last example we present some patterns obtained using the combination of the Picard iteration and the perturbed iteration. We used the Three Ply

Conclusions

In the paper we presented new methods for patterns generation with the use of discrete dynamical systems. The first method uses the perturbation mapping, which alters the orbit of a dynamical system during the iteration process. We also presented some perturbation mappings. The second method uses a combination of the two methods: the standard Picard iteration and the method with the perturbation mapping.

The presented examples show that using the proposed methods we are able to obtain very interesting and diverse patterns. The patterns differ from those obtained with the standard Picard iteration. Comparing to the method from [START_REF] Naud | Automatic Generation of Aesthetic Images for Computer-assisted Virtual Fashion Design[END_REF] the proposed methods give more possibilities of obtaining unrepeatable patterns.

The obtained patterns have an aesthetic value so they can be used as usable patterns, e.g. textile patterns, ceramics patterns, or can be used in jewellery, decoration design or in generating textures using for instance the glyph bombing technique [START_REF] Rost | OpenGL Shading Language, 3rd Edition[END_REF].

  are presented. The parameters for the transformations were following: Hopalong -a = 55, b = -10, c = -42.1, x 0 = 0, y 0 = 0, Zaslavsky -K = 4, q = 7, x 0 = 1, y 0 = 1, Chip -a = 34, b = 1, c = 5, x 0 = 0, y 0 = 0, Quadrup Two -a = 34, b = -1.2, c = -20, x 0 = 0, y 0 = 0, Three Ply -a = 55, b = -10, c = -42.1, x 0 = 0, y 0 = 0.

Fig. 1 .

 1 Fig. 1. The examples of orbits (n = 100000), the top row (from the left): Hopalong, Zaslavsky, Chip, the bottom row (from the left): Quadrup Two, Three Ply

Fig. 2 .

 2 Fig. 2. Orbits of Quadrup Two transformation for different starting points

Fig. 3 .

 3 Fig. 3. Original orbit of Zaslavsky transformation (top) and patterns obtained with perturbation after a fixed number of iterations (bottom)

Fig. 4 .Fig. 5 .

 45 Fig. 4. Original orbit of Chip transformation (top) and patterns obtained with perturbation using different vectors (bottom)

Fig. 6 .

 6 Fig. 6. Original orbit of Three Ply transformation (top-left), pattern obtained with perturbation mapping with a fixed vector (top-right) and patterns obtained using the combination (13) with different values of α (bottom)


Robert Krupiński
email: rkrupinski@wp.pl

Modified Moment Method Estimator for the Shape Parameter of Generalized Gaussian Distribution for a Small Sample Size

Keywords: estimation, generalized Gaussian distribution, moment method

de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The moment method [START_REF] Mallat | A Theory of Multiresolution Signal Decomposition: The Wavelet Representation[END_REF] is very popular method used for the estimation of the power shape parameter of GGD. The most popular approach is when the two moments are set to the values m 1 = 1 and m 2 = 2 (the Mallat's method), but the higher order moments are also used [START_REF] Novey | A complex generalized Gaussian distributioncharacterization, generation, and estimation[END_REF]. Most of the already known estimation methods for the shape parameter of GGD assume that the sample size is large, but in the real case, a set of limited size is only available. A very high value of variance for a very small sample size makes the estimation methods very inaccurate.

The computation of GGD parameters on small sample size is different from on larger ones in a maximum likelihood (ML, [START_REF] Du | Ein sphärisch invariantes Verbunddichtemodell für Bildsignale[END_REF]) framework as was shown in [START_REF] Meignen | On the modeling of small sample distributions with generalized Gaussian density in a maximum likelihood framework[END_REF]. The ML method has a lower variance for a small sample size than the MM method [START_REF] Varanasi | Parametric generalized Gaussian density estimation[END_REF], but it is complex and time consuming. The complexity of ML can be reduced with the One Moment (OM) method [START_REF] Krupiński | Modeling the distribution of DCT coefficients for JPEG reconstruction[END_REF]. The comparison of ML and the Mallat's method for the estimation of GGD shape parameter in the range 0.3 -3 is presented in [START_REF] Krupiński | Approximated fast estimator for the shape parameter of generalized Gaussian distribution[END_REF], where the latter one is not accurate for the whole range of the shape parameter.

The presented approach in this article will not overcome other estimation methods for a small sample size, as it is fully based on the MM framework, but it will enhance the convergence of the MM method for the higher order moments and a small sample size.

In Section 2 the higher order moments are expanded into series in the moment method equation for a small sample size. In Section 3 the numerical results for the integer and real higher order moments are presented.

Material and methods

The probability density function (PDF) of the continuous random variable of GGD [START_REF] Box | Bayesian Inference in Statistical Analysis[END_REF] takes the form

f (x) = λ(p, σ) • p 2 • Γ ( 1 p ) e -[λ(p,σ)•|x-µ|] p , ( 1
)
where Γ (z) = ∫ ∞ 0 t z-1 e -t dt, z > 0, is the Gamma function [START_REF] Olver | Asymptotics and special functions[END_REF], p denotes the shape parameter, µ is the location parameter and λ relates to the variance of the distribution and can be calculated from

λ(p, σ) = 1 σ [ Γ ( 3 p ) Γ ( 1 p ) ] 1 2 , ( 2
)
where σ denotes the standard deviation. The special case of GGD can be observed when the shape parameter equals p = 1 and p = 2, which corresponds to Laplacian and Gaussian distributions respectively. For p → 0, f (x) approaches an impulse function, and for p → ∞, f (x) approaches a uniform distribution.

For µ = 0, the probability density function is centralized around zero.

From the definition of absolute moment

E [|X| m ] = ∫ ∞ -∞ |x| m • f (x)dx, for GGD the following is obtained [5] E m = Γ ( m+1 p ) λ m • Γ ( 1 p ) , (3)
where m can be a real number.

The E m estimated value of the moment can be acquired from the equation

Êm = 1 N • N ∑ i=1 |x i | m (4)
and where N denotes the number of observed variables, and {x 1 , x 2 , . . . , x N } is the collection of N i.i.d zero-mean random variables. Eq. ( 3) for two different moment values m 1 and m 2 and eliminating λ leads to:

Γ ( m2+1 p ) Γ ( 1 p ) m 2 m 1 -1 Γ ( m1+1 p ) m 2 m 1 = E m2 (E m1 ) m 2 m 1 , ( 5
)
which is the method for the estimation of the shape parameter based on two moments.

The moments (Eq. ( 4)) in the rhs of Eq. ( 5) can be expanded into series. Each component |x i | m in the sum (Eq. ( 4)) can be represented in the form (1 ± y) m , which can be expanded into series for m > 0 and |y| ≤ 1. For x i ∈< 0, 1 >, it can be written (by Binomial Theorem)

(x i ) m = (1 + (x i -1)) m = 1 + m • (x i -1) + m•(m-1) 2! • (x i -1) 2 + m•(m-1)•(m-2) 3! • (x i -1) 3 + . . . . ( 6
)
For

x i ∈< -1, 0), it takes the form (-x i ) m = (1 -(x i + 1)) m = 1 -m • (x i + 1) + m•(m-1) 2! • (x i + 1) 2 - m•(m-1)•(m-2) 3! • (x i + 1) 3 + . . . . ( 7
)
For x i > 1, it can be written

(x i ) m = |⌊x i ⌋| m • (1 + ( xi |⌊xi⌋| -1 ) ) m = |⌊x i ⌋| m • { 1 + m • ( xi |⌊xi⌋| -1 ) + m•(m-1) 2! • ( xi |⌊xi⌋| -1 ) 2 + m•(m-1)•(m-2) 3! • ( xi |⌊xi⌋| -1 ) 3 + . . . } . ( 8
)
The operator ⌊x⌋ rounds the elements of x toward zero, resulting in an integer. For x i < -1, it can be written

(-x i ) m = |⌊x i ⌋| m • (1 - ( xi |⌊xi⌋| + 1
) ) m = |⌊x i ⌋| m • { 1 -m • ( xi |⌊xi⌋| + 1 ) + m•(m-1) 2! • ( xi |⌊xi⌋| + 1 ) 2 - m•(m-1)•(m-2) 3! • ( xi |⌊xi⌋| + 1 ) 3 + . . . } . ( 9
)
The equations ( 6)-( 9) are used in the approximation of the moment (Eq. ( 4)) and used in Eq. ( 5) for the solution of p. Taking into account only a limited number of series components, the influence of outliers on the MM equation can be alleviated for a small sample size and the higher order moments and consequently reduce RMSE and improve the convergence.

Calculation

The equations from the article were validated with the GGD generator with fixed variance to unity and varying shape parameter in the range p ∈< 0.3, 6 >.

The relative mean square error was applied for the comparison of the estimators output. RMSE was calculated from the equation

RM SE = 1 M M ∑ i=1 (p -p) 2 p 2 , ( 10
)
where p is a value estimated by the model and p is a real value of a shape parameter. M denotes the number of repetitions.

A small sample size may lead to difficulties with the root finding of the MM method (Eq. ( 5)), therefore, the stop condition was set to tolX = 1e -4 and tolY = 1e -4. When at least one solution is not available for the specified conditions, then the plots have the missing points for such a case. Fig. 1 depicts RMSE for the MM method of the estimation of the shape parameter p of GGD with a varying shape value p and a selected fixed sample size N = 10 4 . Each iteration was repeated M = 10 3 . For a relatively large sample size, RMSE is not constant for all shape parameters of GGD for the MM method. The RMSE value depends on the selected moments m 1 and m 2 and is relatively small. As an example, the integer moments were selected m 1 = 4, m 2 = 8 and the real ones: m 1 = 3.5, m 2 = 5.5. The higher order moments exhibit the better RMSE performance of MM for the higher values of p. In the case of a small sample size, the increase of RMSE as well as the reduction of a convergence range are observed (Fig. 2). The higher order moments assure convergence of the MM equation for the range with higher values of power shape parameter of GGD, whereas the Mallat's method may not have a solution there. The range of convergence of the MM equation for the Mallat's method also reduced. Fig. 2 depicts the plots for N = 100 and N = 200 and each iteration was repeated M = 10 4 . The same repetition count is kept for the rest of calculations. By the application of Eq. ( 6)-( 9) the influence of outliers on the MM equation can be alleviated for a small sample size and the higher order moments and the range of convergence can be extended toward the smaller sample sizes. The expansion of |x i | m in Eq. ( 4) up to and including a k-th component, i.e., ( xi ⌊xi⌋ ± 1) k will be denoted as E m = f (k). Fig. 3 depicts that selecting only a few components in the moments calculation (Eq. ( 6)-( 9)) RMSE is stable and smaller than in the usual calculation. The selected expansions [START_REF] Box | Bayesian Inference in Statistical Analysis[END_REF] are plotted for the integer moments m 1 = 4, m 2 = 8 and for the shape parameter p = 6 and exhibit better RMSE performance below N = 300.

E m1 = f (3), E m2 = f (6) and E m1 = f (3), E m2 = f
The same behavior for the integer moments m 1 = 4 and m 2 = 8 and for the shape parameter p = 5 is observed. For the selected expansions E m1 = f (3),

E m2 = f (6) and E m1 = f (3), E m2 = f (8)
, RMSE is stable and smaller than in the usual calculation below N = 260 (Fig. 4).

The selected expansions E m1 = f (3), E m2 = f (6) and E m1 = f (3), E m2 = f (8) and the MM method for the integer moments m 1 = 4 and m 2 = 8 for a small sample size N = 169 and for a wider range of p are compared in Fig. 5. It can be noticed that the expansions E m1 = f (3), E m2 = f (6) and E m1 = f (3), E m2 = f (8) resulted in smaller RMSE and assured convergence for wider range of p.

In the case of a small sample size for the real higher order moments, the convergence of the MM equation can be also improved toward the observation smaller in size. Fig. 6 depicts a few components in the moments calculation (Eq. ( 6)-( 9)) comparing to the whole moment calculation (m 1 = 3.5, m 2 = 5.5) over a range of varying sample size N for the fixed value p = 6 in the GGD generator. It can be also noted that RMSE is longer stable for E m1 = f (3) and E m2 = f (6) toward smaller N and has a smaller value too.

For the moments m 1 = 3.5 and m 2 = 5.5 and for the shape parameter p = 5, one of the expansions E m1 = f (3) and E m2 = f (6) has smaller RMSE and is stable for longer range of N than the full moments, whereas another expansion E m1 = f (3) and E m2 = f (4) gains RMSE smaller and is still stable below N = 250 (Fig. 7).

The selected expansions Comparison of RMSE for the MM method (m1 = 3.5, m2 = 5.5) and selected expansions of the estimation of the shape parameter p of GGD with a varying sample size N and for the selected fixed value p = 5 in the GGD generator small sample size N = 300 and for a wider range of p are compared in Fig. 8. It can be noticed that the expansions resulted in smaller RMSE and assured convergence for wider range of p. For the one of the expansions, i.e., E m1 = f (3) and E m2 = f (4) the increase of RMSE is observed comparing to the normal moments calculation, but the results are stable for the wider p range. The same approach can be applied for other higher order moments both integer and real, but not all expansions lead to better performance.

E m1 = f (3), E m2 = f (4) and E m1 = f (3), E m2 = f (6)

Conclusions

The article focuses on the estimation of the shape parameter of GGD when a small size is only available. The experimental results and the comparison with already known MM method, which show that the method proposed is very useful from a practical point of view. The smaller sample size the variance of the MM estimator increases and the moment method equation may not converge to a real solution for some sample sets. For the higher order moments, the presented approach allows to lower RMSE and ensure a real solution for the smaller sample sets of GGD for the modified MM method.

Fig. 1 .

 1 Fig. 1.Comparison of RMSE for the MM method of the estimation of the shape parameter p of GGD with a varying shape value p and a selected fixed sample size N = 10 4 (a relatively large sample size) for selected moments

Fig. 2 .

 2 Fig. 2. Comparison of RMSE for the MM method of the estimation of the shape parameter p of GGD with a varying shape value p and a selected fixed sample size N = 100 and N = 200 for selected moments

Fig. 3 .

 3 Fig. 3.Comparison of RMSE for the MM method (m1 = 4, m2 = 8) and selected expansions of the estimation of the shape parameter p of GGD with a varying sample size N and for the selected fixed value p = 6 in the GGD generator

Fig. 4 .Fig. 5 .Fig. 6 .

 456 Fig. 4.Comparison of RMSE for the MM method (m1 = 4, m2 = 8) and selected expansions of the estimation of the shape parameter p of GGD with a varying sample size N and for the selected fixed value p = 5 in the GGD generator

Fig. 8 .

 8 Fig.8. Comparison of RMSE for the MM method (m1 = 3.5, m2 = 5.5) and selected expansions of the estimation of the shape parameter p of GGD with a varying shape value p and a selected fixed sample size N = 300


Rahma Boucetta
email: rboucetta@yahoo.fr

Generalized Predictive Control for a Flexible Single-link Manipulator

Keywords: Flexible manipulator robot, dynamic model, generalized predictive control, fuzzy supervisory control

This paper presents the development of a generalized predictive controller applied to a flexible single-link manipulator robot to compare to a fuzzy supervisory controller in input tracking and end-point vibration suppression. A dynamic model of the flexible manipulator is derived using finite elements method and Lagrange's equations to determine dynamics behavior. A generalized predictive controller is then developed and introduced in the system closed-loop to minimize end-point residual vibrations. A fuzzy supervisory controller is also synthesized to compare simulation results between the two methods of control in terms of input tracking and disturbance rejection.

Introduction

Robotic manipulators are generally built using heavy material to maximize stiffness, in an attempt to minimize system vibration and achieve good positional accuracy. As a consequence, such robots are usually heavy with respect to the operating payload. The operation speed of the robot manipulation is limited, so the actuators size is increased boosting energy consumption and increasing the overall cost. Moreover, the robot has a low payload to robot weight ratio. In order to solve these problems, robotic manipulators are designed to be lightweight. Conversely, flexible manipulators exhibit many advantages over their rigid counterparts: they require less material, are lighter in weight, have higher manipulation speed, lower power consumption, require smaller actuators, are more maneuverable and transportable, are safer to operate due to reduced inertia, have enhanced back-drive ability due to elimination of gearing, have less overall cost and higher payload to robot weight ratio. However, the control of flexible robotic manipulators to maintain accurate positioning is an extremely challenging problem. Due to the flexible nature and distributed characteristics of the system, the dynamics are highly non-linear and complex. Problems arise due to precise positioning requirement, vibration due to system flexibility, difficulty in obtaining an accurate model and non-minimum phase characteristics of the system. Therefore, flexible manipulators have not favored in production industries, due to un-attained end-point positional accuracy requirements in response to input commands. Thus, the design of control algorithms for flexible systems possessing nonlinear time-varying and ill-modeled dynamics presents great challenges for all conventional methodologies. M.A. Arteaga and B. Siciliano collected a number of recent results on modeling, nonlinear control and observer for flexible-link manipulators [START_REF] Wang | Advanced Studies of Flexible Robotic Manipulators, Modeling, Design, Control and Applications[END_REF]. S.S. Ge proposed energy-based robust control strategies for the control of flexible link robots without using the dynamics of the systems explicitly [START_REF] Wang | Advanced Studies of Flexible Robotic Manipulators, Modeling, Design, Control and Applications[END_REF]. O. Al Jarrah, Y.F. Zheng and K.-Y.

Yi presented three approximation methods of the optimal trajectories and a compliant control scheme [START_REF] Wang | Advanced Studies of Flexible Robotic Manipulators, Modeling, Design, Control and Applications[END_REF]. R.N. Banavar and P. Dominic applied an LQG/H  controller for a flexible manipulator [START_REF] Banavar | An LQG/H Controller for a Flexible Manipulator[END_REF]. S.B. Choi and J.W. Cheon proposed a vibration control of a single-link flexible arm subjected to disturbances [START_REF] Choi | Vibration control of a single-link flexible arm subjected to disturbances[END_REF]. A neural network control is developed by C.-F.J. Kuo and C.-J. Lee for a rotating elastic manipulator [START_REF] Kuo | Neural Network Control of a Rotating Elastic Manipulator[END_REF]. L. Tian and C. Collins proposed firstly a dynamic recurrent neural network-based controller for a rigid-flexible manipulator system [START_REF] Tian | A dynamic recurrent neural network-based controller for a rigidflexible manipulator system[END_REF] and secondly an adaptive neuro-fuzzy control for a flexible manipulator [START_REF] Tian | Adaptive neuro-fuzzy control of a flexible manipulator[END_REF]. A self-organizing fuzzy logic controller is used by G.L.C.M. de Abreu and J.F. Ribeiro for the active control of flexible structures using piezoelectric actuators [START_REF] De Abreu | A self-organizing fuzzy logic controller for the active control of flexible structures using piezoelectric actuators[END_REF]. Soft computing methods are applied to the control of a flexible robot manipulator by B. Subudhi and A.S. Morris [START_REF] Subudhi | Soft computing methods applied to the control of a flexible robot manipulator[END_REF]. It is in the last area that our contribution can be located.

The main purpose of this work is to determine the adequate controller to minimize human intervention and to increase performance responses. So, this paper is organized as follows: Section 2 describes the flexible manipulator system and how to derive its dynamic model using Lagrange's equations and finite elements method. The openloop analysis is given in section 3. Section 4 presents the GPC controller applied to the flexible single-link manipulator. Section 5 developed the supervisory fuzzy controller given to the system. Finally, a comparative assessment of performances between the different strategies in terms of vibration suppression and input tracking is presented and discussed.

The Flexible Manipulator System

A schematic representation of the single-link flexible manipulator system is shown in figure 1, where a control torque τ(t) is applied at the hub of a motor with E, I, ρ, L and I H represent Young's modulus, second moment of area, mass density per unit volume, length, and hub inertia moment respectively [START_REF] Meirovitch | Elements of Vibration Analysis[END_REF].

The angular displacement of the link in the X 0 OY 0 coordinates is denoted by θ(t). w(x,t) represents the elastic deflection of the manipulator at a distance x from the hub, measured along the OX axis. X 0 OY 0 and XOY represent the stationary and moving frames respectively.

Fig. 1. Flexible manipulator scheme

The height (width) of the link is assumed to be much greater than its depth, thus allowing the manipulator to vibrate dominantly in the horizontal direction (X 0 OY 0 plane). To avoid difficulties arising from time varying lengths, the length of the manipulator is assumed to be constant. Moreover, the shear deformation, the rotary inertia and the effect of axial force are ignored. For an angular displacement θ and an elastic deflection w, the total displacement y(x,t) of a point along the manipulator at a distance x from the hub can be described as a function of both the rigid body motion θ(t) and the elastic deflection w(x,t), i.e. [3, 7]

( ) ( ) ( ) (1)
Thus, by allowing the manipulator to be dominantly flexible in the horizontal direction, the elastic deflection of the manipulator can be assumed to be confined to the horizontal plane only. Kinetic energy of the flexible manipulator, depending of hub rotation, modes rotation in the X0OY0 and XOY frames, has the following expression

̇ ( ̇ ̇) ( ̇ ̇) (2)
Potential energy just depending of link flexibility has the form

After applying Lagrange's equations, the dynamic model can be written as

( ) ̈ ̈ (4) ̈ ̈ (5)
where M, K and L are the mass matrix, the stiffness matrix and the length array respectively, and q is the elastic modes vector.

Dynamic Behavior

The dynamic equations can be presented in a state-space form as

̇
where the state-space matrices are

( ) ( ) ( ) ( )
The state and control vectors are given by

( ̇ ̇ ̇ ) ( )
In order to simulate the flexible manipulator system, an aluminum link of dimensions L=0.61m and S=3×10 -5 m², with E=200×10 9 N/m², I=2.5×10-12 m 4 , I H =4.3×10 -3 Kg.m² and ρ=7.8×103Kg/m 3 is considered [START_REF] Boucetta | Building and Control of a Flexible One-link Manipulator[END_REF]. The link is discretized into two elements. Solving the state-space matrices gives the vector of states v, that is, the hub angle, the elastic modes and their velocities. The derived dynamic model is a nonminimum phase system, not strictly proper, and unstable. Also, the model has zeros very close to the imaginary axis; this deteriorates the time domain performance of the closedloop system. Generally, linear models of flexible structures used in design of controllers are derived under restrictive assumptions which are often not valid for large motions that occur during slewing maneuvers. Hence, considerable uncertainty in the linear model exists. Another feature characteristic of lightly damped systems is the occurrence of poles (±159.67j, ±37j, 0, 0) and zeros (±158j, ±25j) very close to the imaginary axis that gives rise to ill-conditioned systems. The state-space matrices arising out of such systems have largely separated singular values, posing considerable computational difficulty in controller design. In the spectral density given by figure 2, the vibration frequencies of the system are obtained as 37rad/s and 160rad/s, i.e. 5.9Hz and 25.46Hz, and the magnitude of frequency response for the two resonance modes are 122dB and 68.8dB.

The prediction margin ranged from to , and the system delay is known a priori, so, the ( ) first lines of matrices are not taken into account, and eq. 9 and eq. 10 are rewritten as

̂ ̂

The aim of control is to compute optimal control sequence minimizing the cost function. Length of this vector is reduced to ( ), and ( ) last columns of matrix are not taken into account. Consequently, a new matrix and a new expression of ̂ ̂ ̂ ̂ ̂

The cost function has the following matrix form

̂ [ ̂ ]
The minimal value of the cost function is

[ ̂ ]
The optimal control for the minimal cost function has the following form with [ ] Finally, the resulted control signal is the first element of the control vector , given by the following expression

( ) ( ) ( ) with
( )
Fig. 3. GPC control scheme

Fuzzy Supervisory Control

A highest level supervisor uses any available data from the control system to characterize the system's current behavior so that it knows how to change the controller and ultimately achieve the desired specifications. In addition, the supervisor can be used to integrate other information into the control decision-making process. It can incorporate certain user inputs, or inputs from other subsystems [START_REF] Passino | Fuzzy Control[END_REF]. Conceptually, the design of the supervisory controller can then proceed in the same manner as it did for direct fuzzy controllers [START_REF] Abdelkrim | Compensation des mouvements respiratoires par la commande prédictive généralisée à modèle répétitif d'un robot médical[END_REF]: either via the gathering of heuristic control knowledge or via training data that we gather from an experiment. The type of heuristic knowledge that is used in a supervisor may take one of the following two forms  Information from a human control system operator who observes the behavior of an existing control system and knows how this controller should be tuned under various operating conditions.  Information gathered by a control engineer who knows that under different operating conditions controller parameters should be tuned according to certain rules.

A higher level of control can be achieved for monitoring and adjusting the direct fuzzy controller. The expert controller expands or compresses the universes of discourse by simply changing the scaling gains. When the universe is expanded, a coarse control given by the table 1 is achieved, and when it is compressed, a fine control given by the table 2 is achieved. The supervisor would be a fuzzy system that can gradually rather than abruptly switch between the two conditions using a Sugeno fuzzy system based on the following rule-base:

 If error is negative (positive) Then control action is coarse  If error is zero Then control action is fine Table 1. Rule-base for the coarse control u ε' -3 -2 -1 0 +1 +2 +3 ε -3 -4 -3 -3 -2 -2 -1 0 -2 -3 -3 -2 -2 -1 0 +1 -1 -3 -2 -2 -1 0 +1 +2 0 -2 -2 -1 0 +1 +2 +2 +1 -2 -1 0 +1 +2 +2 +3 +2 -1 0 +1 +2 +2 +3 +3 +3 0 +1 +2 +2 +3 +3 +4
u ε' -3 -2 -1 0 +1 +2 +3 ε -3 -4 -4 -3 -3 -2 -1 0 -2 -4 -3 -3 -2 -1 0 +1 -1 -3 -3 -2 -1 0 +1 +2 0 -2 -1 0 0 0 +1 +2 +1 -2 -1 0 +1 +2 +3 +3 +2 -1 0 +1 +2 +3 +3 +4 +3 0 +1 +2 +3 +3 +4 +4

Simulation Results

To study the GPC and supervisory control performances and to compare simulation results in terms of input tracking, vibration suppression and disturbance rejection, the GPC controller is introduced in the closed-loop position control of the flexible singlelink manipulator. Sampling time is chosen 0.01s. Figure 5 depicts the time evolutions of both consign and hub angle responses. A transient appears at t=10s, corresponding to the start of the step consign chosen . The control signal can compensate the fast change of consign at t=10s. So, the hub angle output varies before the consign signal and has very fast rise explaining the appearance of an overshoot peak over the consign value. After that, the plant output oscillates around the desired value with very low vibration magnitude.

Fig. 4. Consign and output responses without disturbances

In figure 6, an impulse disturbance with amplitude of is applied at t=15s, the plant output receives a small variation and continues with the same specifications. 7 gives the hub angle response deduced from the fuzzy supervisory control, the output response has a fast transient phase to reach consign without overshoot, with rise time and settling time about 0.5s. The introduction of a disturbance at t=15s in the fuzzy supervisor control case yields the result of figure 8, we can notice that the hub angle receives a slight peak and returns to its stable state. For the two controllers, we can easily noticed that the GPC controller gives a fast hub angle response with a good input tracking but it has an overshoot in the rise, contrary to the fuzzy supervisory controller that has a slower rise time compared to the GPC, without vibrations or overshoot. The GPC controller has a good disturbance rejection against the fuzzy supervisor. A comparative assessment between the GPC and the fuzzy supervisory control strategies applied to a flexible single-link manipulator has been presented in this paper. A dynamic model of the flexible manipulator system is first derived using Lagrange's equations and finite element method. A GPC controller is then introduced in the closed-loop of the flexible system. Next, a fuzzy supervisory controller is developed to be incorporated in the system closed loop to increase level of desired performances. Simulation results are compared in terms of input tracking, disturbance rejection and vibration reduction.

Fig. 5 .Fig. 6 .Fig. 7 .

 567 Fig. 5. Consign and output responses with disturbance at t=15s

Table 2 .

 2 Rule-base for the fine control

The Generalized Predictive Controller was introduced by Clarke and al. in 1987 [START_REF] Clarke | Properties of Generalized Predictive Control[END_REF].

The computation of the output predictions supposes the knowledge of a model of the system that is in the ARIMAX form:

where:

 is the backward operator  ( ) is the output signal  ( ) is the input signal  ( ) is the disturbance (white noise) process with ( ( ))  A and C are monic polynomial matrices where C can be used to model a colored noise. To simplify the problem, we consider that ( ) where c is a monic polynomial  B is a polynomial matrix  The operator is defined as used to make noise be non-stationary, which suitable to model any perturbation in a control loop.

The polynomial matrices A, B, C are respectively of order , , and :

where d is a positive integer representing the delay of the system. The operator ( ) allows the rejection of constant perturbations and is equivalent to the introduction of an integral action in the controller.

Cost Function

The GPC is based on the minimization of a cost function J over a finite receding horizon [START_REF] Abdelkrim | Compensation des mouvements respiratoires par la commande prédictive généralisée à modèle répétitif d'un robot médical[END_REF][START_REF] Abdelkrim | commande par retour visuel d'un robot médical[END_REF][START_REF] Ginhoux | Compensation des mouvements physiologiques en chirurgie robotisée par commande prédictive[END_REF]:

with and

where λ is a positive scalar and , and are positive integers defined as follows:

 is the minimum costing horizon,  is the maximum costing horizon,  is the length of the control cost horizon,  weights the relative importance of the control energy,  is the sampling time  ( ) is the reference trajectory The aim is to compute the future control increments, so as to drive the actual system outputs towards the theoretical ones, or, equivalently, to compensate for the measurement disturbances. This controller is predictive because it takes into account the future references. Indeed, the minimization of (6) requires the computation of predictions of the output using the future reference signals. The arguments of the minimization are the future steps of the control input.

Control computing

The minimization of the cost function on a finite horizon allows us to compute the predictor on j optimal number of sampling periods to predict real output in ( ) sample periods. The vector is setting as

The prediction depending of past and present measures ̂ is defined as

Similarly, the prediction depending of the sequence of future control increments ̂ minimizing the cost function is given by ̂ [START_REF] Subudhi | Soft computing methods applied to the control of a flexible robot manipulator[END_REF]


Kumar Sanjoy

Ruhul Paul

Daryl Sarker
email: r.sarker@adfa.edu.au

Essam
email: d.essam@adfa.edu.au

A Disruption Recovery Model in a Production-Inventory System with Demand Uncertainty and Process Reliability

Keywords: Production inventory, series of disruptions, demand uncertainty, process reliability, genetic algorithm

This paper develops a risk management tool for a productioninventory system that involves an imperfect production process and faces production disruption and demand uncertainty. In this paper, the demand uncertainty is represented as fuzzy variable and the imperfectness is expressed as process reliability. To deal with the production scheduling in this environment, a nonlinear constrained optimization model has been formulated with an objective of maximizing the graded mean integration value (GMIV) of the total expected profit. The model is applied to solve the production-inventory problem with single as well as multiple disruptions on a real time basis that basically revises the production quantity in each cycle in the recovery time window. We propose a genetic algorithm (GA) based heuristic to solve the model and obtain an optimal recovery plan. A numerical example is presented to explain usefulness of the developed model.

Introduction

Batch production is a well accepted technique in advanced manufacturing and logistics management system. Production lot size is determined to minimize the costs of the system. There are numerous industries, such as pharmaceutical, textile and food, that produce the products using the batch production systems. There are several risks factors in real life problems which should be taken into consideration when production system is analyzed. Production disruptions i.e. raw material shortage, machine breakdown, labor strike, or any other production interruptions are very common scenario in the production systems. Moreover, it is very difficult to find the production process that produces 100% non-defective products. So the production process reliability, can be less than 100%, is also an important factor because of the imperfect production process. In real life situations, product demand cannot be known with uncertainty. In this paper, the process reliability and demand uncertainty are consid-ered with production disruption to make the research problem very close to the practical scenario. Over the last few decades, one of the most widely studied research topics, in operations research and industrial engineering, is the production inventory system. Few examples of such studies in single stage production inventory system include a singleitem inventory system with non-stationary demand process [START_REF] Graves | A single-item inventory model for a nonstationary demand process[END_REF], determination of lot size and order level for a single-item inventory model with a deterministic timedependent demand [START_REF] Dave | A deterministic lot-size inventory model with shortages and a linear trend in demand[END_REF], a single-item periodic review stochastic inventory system [START_REF] Chan | A dynamic analysis of the single-item periodic stochastic inventory system with order capacity[END_REF] and a single-item single-stage inventory system with stochastic demand in a periodic review where the system must order either none or at least as much as a minimum order quantity [START_REF] Kiesmüller | Single item inventory control under periodic review and a minimum order quantity[END_REF].

The above studies, with many others, are conducted under ideal conditions. However, production disruption is a very familiar event in the production environment. Production disruption is defined as any form of interruption that may cause due to shortage of material, machine breakdown and unavailability, or any other form of disturbance. The development of an appropriate recovery policy can help to minimize the loss and maintain the goodwill of the company. Lin and Gong [START_REF] Lin | On a production-inventory system of deteriorating items subject to random machine breakdowns with a fixed repair time[END_REF] analysed the impact of machine breakdown on EPQ model for deteriorating items in a single stage production system with fixed period of repair time. Widyadana and Wee [START_REF] Widyadana | Optimal deteriorating items production inventory models with random machine breakdown and stochastic repair time[END_REF] extended the model of Lin and Gong [START_REF] Lin | On a production-inventory system of deteriorating items subject to random machine breakdowns with a fixed repair time[END_REF] for deteriorating items with random machine breakdown and stochastic repair time with uniform and exponential distribution. A disruption recovery model for single stage and single item production system is developed by Hishamuddin et al. [START_REF] Hishamuddin | A disruption recovery model for a single stage production-inventory system[END_REF] and the model was formulated for a single disruption, for recovering within a given time window, considering back order as well as lost sales option. Recently, a transportation disruption recovery model in a two-stage production and inventory system was developed by Hishamuddin et al. [START_REF] Hishamuddin | A Recovery Model for a Two-Echelon Serial Supply Chain with Consideration of Transportation Disruption[END_REF]. In the production and inventory modelling, numerous studies have been performed considering supply disruptions. Parlar and Perry [START_REF] Parlar | Inventory models of future supply uncertainty with single and multiple suppliers[END_REF] developed inventory models considering supplier availability with deterministic product demand under the continuous review framework. Özekici and Parlar [START_REF] Özekici | Inventory models with unreliable suppliers in a random environment[END_REF] considered back orders to analyse a productioninventory model under random supply disruptions modelled as a Markov chain. Recently, other models of supply disruptions considering deterministic product demand in the inventory models have been studied [START_REF] Mohebbi | An inventory model with non-resuming randomly interruptible lead time[END_REF] and [START_REF] Qi | A continuous-review inventory model with disruptions at both supplier and retailer[END_REF].

There are some recent studies, where reliability of the imperfect production process has been considered. At first, process reliability is considered by Cheng [START_REF] Cheng | An economic production quantity model with flexibility and reliability considerations[END_REF] in a single period inventory system and formulated as unconstrained geometric programming problem. Later, it was extended in [START_REF] Bag | A production inventory model with fuzzy random demand and with flexibility and reliability considerations[END_REF] by considering fuzzy random demand. Later, process reliability of the imperfect production process was incorporated to determine the optimal product reliability and production rate that achieved the biggest total integrated profit [START_REF] Sana | A production-inventory model in an imperfect production process[END_REF], to study unreliable supplier in a single-item stochastic inventory system [START_REF] Mohebbi | An inventory model with non-resuming randomly interruptible lead time[END_REF] and to analyze an EPQ model with price and advertising demand pattern under the effect of inflation [START_REF] Sarkar | An inventory model with reliability in an imperfect production process[END_REF].

Many authors considered only fuzzy characteristics of the variables to tackle uncertainty. Lee and Yao [START_REF] Lee | Economic production quantity for fuzzy demand quantity, and fuzzy production quantity[END_REF] introduced fuzzy senses in the EPQ model considering demand and production per day as fuzzy variables. Later fuzzy product quantity as a triangular fuzzy number [START_REF] Chang | Fuzzy production inventory for fuzzy product quantity with triangular fuzzy number[END_REF], order quantity and total demand quantity as triangular fuzzy numbers Yao et al. [START_REF] Yao | Fuzzy inventory without backorder for fuzzy order quantity and fuzzy total demand quantity[END_REF], demand as a fuzzy random variable in a single-period inventory model [START_REF] Dutta | A single-period inventory model with fuzzy random variable demand[END_REF] are considered in developing production inventory modelling. Recently, Islam and Roy [START_REF] Islam | Fuzzy multi-item economic production quantity model under space constraint: A geometric programming approach[END_REF] developed a modified geometric programming program in an EPQ model under storage space constraint and reliability of the production process considering inventory related costs, storage spaces and others parameters as triangular fuzzy number.

In the previous studies of production inventory modelling, no study considered demand uncertainty and process reliability to develop a disruption recovery model. Also most of the previous studies focused on developing recovery plan only from single disruption. In this paper, a real time disruption recovery model is developed where the production process faces single or multiple disruptions. Other risk factors, process reliability and demand uncertainty, are also incorporated to make the model realistic. Finally, a genetic algorithm based heuristic is proposed to solve the model with single or multiple disruptions on a real time basis.

Problem Definitions

In the real life production system, disruptions are very common scenario and it can happen at any time at any point. It needs to develop an optimal plan to recover from those production disruptions. Revision in the production quantities and use of the idle timeslot of the systems are the significant ways to obtain the recovery plan [START_REF] Hishamuddin | A disruption recovery model for a single stage production-inventory system[END_REF]. After each disruption, production quantities in each cycle during the recovery period are revised. We develop a solution approach to obtain the recovery plan that deals with single or series of disruptions on a real time basis.

Fig. 1. Disruption recovery plan

The recovery plan after a production disruption is presented in Figure 1. The recovery plan is a new schedule which includes the revised production quantities in each cycle to maximize the total profit in the recovery period. Number of cycles allocated to return to the original production schedule from the disrupted cycle is known as recovery period. The first disrupted production cycle is considered as first cycle (l=1). Now an optimal recovery plan is proposed to revise the production quantities ( ) to recover from that disruption which is shown as black dashed line in the Figure 1. Again after the second disruption within the recovery period of previous disruption, production quantities in each cycle ( ) are revised by considering the effect of both disruptions which is shown as red dashed line In this study, we have made a number of assumptions as follows.

i. Production rate is greater than GMIV of the demand rate. ii.

Single item is produced in the system. iii.

All products are inspected and defective products are rejected. iv.

Total cost of interest and depreciation per production cycle ( ) is inversely related to set-up cost ( ) and is directly related to the process reliability (r) according to following general power function [START_REF] Cheng | An economic production quantity model with flexibility and reliability considerations[END_REF]:

( ) Where a, b and c are positive constants chosen to provide the best fit of the estimated cost function.

Notations used in the Study

The In this section, economic production lot size ( ), equations for different costs and revenues and final objective function are derived for the single stage production inventory system that considers process reliability and demand uncertainty. Economic lot size is calculated to minimize the total annual set-up and holding cost. For a single item production system, with lot-for-lot condition under ideal situation [START_REF] Sarker | An optimal batch size for a production system operating under periodic delivery policy[END_REF] with process reliability, the economic production quantity can be formulated as:

√ (1)

Costs and Revenues Formulation

Holding, set-up, back order, lost sales, production, rejection, inspection and depreciation costs are identified as the relevant costs. Holding cost is determined as unit holding cost multiplied by total inventory during the recovery period which is equivalent to the area under the curve of Figure 1. Set-up cost is calculated as cost pet set-up multiplied by number of set-up in the recovery period. Back order cost is determined as unit back order cost multiplied by back order units and it's time delay [START_REF] Hishamuddin | A disruption recovery model for a single stage production-inventory system[END_REF]. Lost sales cost is determined as unit lost sales cost multiplied by lost sales units [START_REF] Hishamuddin | A disruption recovery model for a single stage production-inventory system[END_REF]. Unit production cost multiplied by total quantity produced during the recovery period is the total production cost. Rejection cost is determined as unit rejection cost multiplied by total rejected quantities [START_REF] Paul | Development of a production inventory model with uncertainty and reliability considerations[END_REF]. Inspection cost is considered as a certain percentage of the production cost [START_REF] Paul | Development of a production inventory model with uncertainty and reliability considerations[END_REF]. Cost of interest and depreciation equation is considered as a general power function [START_REF] Cheng | An economic production quantity model with flexibility and reliability considerations[END_REF]. The model is generalized by considering the production quantity in cycle i after the n th disruption as and the original production quantity as .

Holding cost

[ ( ) ( ) ∑ ( ) ] (2) Set-up cost (3) Production cost (∑ ) (∑ ) (4) Rejection cost ( ) (∑ ) ( ) (∑ ) (5) Inspection cost (∑ ) (6) Cost of interest and depreciation ( ) ( ) (7) Back-order cost [( ) [ ] ∑ [ ( ) ∑ ∑ ( ) ]] (8) Lost sales cost (∑ ∑ ) (9
) Selling price of the acceptable items, which is revenue, in the recovery period is determined as unit selling price multiplied by the demand in the recovery period [START_REF] Paul | Development of a production inventory model with uncertainty and reliability considerations[END_REF].

Revenues ̃[∑ ] (10)
In this paper, demand data are considered as a triangular fuzzy number (TFN). Demand TFN and associated probabilities are taken as a triplet ( ̅̅̅ ) and ( ̅̅̅) respectively. Where, k= 1, 2, 3...., v. Then the fuzzy expected profit function will also be a TFN, ̃ ( ̅̅̅̅ ) which is determined as follows:

( )] ∑ ( ) ] ∑ ( )] ∑ ̅̅̅̅ ̅̅̅̅
Here the -level set of the fuzzy number ̃ are considered as ( )

( )] [ ( ( )) ( ( ))]
and different -cut intervals for the fuzzy number ̃ are obtained for different between 0 and 1.Taking, -cut on both sides of equation of ̃. ̃ ∑ ̃ ̃ The arithmetic interval of fuzzy demand and associated probabilities using an -cut is determined as follows.

̃ [ (

By using these arithmetic intervals, ̃ is evaluated as:

̃ [[ ∑ [ ( )] [ ( )] ] [ ∑[ ( )] ( )] ]]
From the representation of graded mean integration methods based on the integral value of the graded mean -level of the LR-fuzzy number of the total profit, ( ) and ( ) are obtained as follows.

(

] The unique fuzzy number ( ̃) is determined by substituting the value of ( ) and ( ) to the equation ( 13),

( ̃) ∫ [ [ ∑ [ ( )] [ ( )] ]] ∫ [ [ ∑[ ( )] ( )] ]]
Population size: 100; Population type: Double vector; Crossover fraction: 0.8; Maximum number of generations: 3000; Function tolerance: 1e-6; Non linear constraint tolerance: 1e-6 and other parameters are set as default of the optimization toolbox.

Results Analysis

Results have been analysed for both single and multiple disruptions on a real time basis. For single disruption, there is only one random disruption in the system and there is no more disruption within the recovery period. For multiple disruptions, there is a series of disruptions, one after another, on a real time basis in the system.

Results Analysis for Single Disruption

Following data are considered to analyze the results for single disruption:

, and demand data are considered as TFN which is shown in the Table

Results Analysis for a Series of Disruptions

In this case, a series of disruptions on a real time basis is considered, which is shown in Table 3. In this series of disruptions, seven dependent disruptions are considered and each one occurs within the recovery period of the previous disruption. Other data remain same as in section 7.1. 3 0 0.0078 The production system with multiple disruptions is also solved using the GA based heuristic on a real time basis. The results are obtained from 30 different runs. Production quantity in each cycle is revised after each disruption considering the effect of entire dependent disruptions to maximize the total profit in the recovery period. The best recovery plan obtained from the heuristic for the series of disruptions is shown in Table 4. The production system returns to original schedule from the sixth cycle after each disruption with , and so on.

Conclusions

The objective of this research was to incorporate demand uncertainty and process reliability in developing a disruption recovery model for managing risk in a production inventory system. A single or a series of disruptions on a real time basis was considered to make the model applicable in practical problems. The model was formulated as a non-linear constrained optimization problem and generalized by formulating the model for the n th disruption. A genetic algorithm based heuristic was proposed to solve the model with single or multiple disruptions on a real time basis. This model can be applied in an imperfect production process where the process countenances a single or multiple production disruptions and product demand is uncertain. The model can be extended by considering multiple stages in the production system.

1 .

 1 It will be continued same way if there is any other disruption. The model is generalized by formulating for the n th disruption.

Fig. 2 .

 2 Fig. 2. Flowchart of proposed GA based heuristic

  Productions up time for cycle i in the recovery period after n th disruption Unit back order cost per unit time ($ per unit per unit time)

	Model Formulation
	Unit lost sales cost ($ per unit)
	Per unit production cost ($ per unit)
	Rejection cost per unit ($ per unit)
	Inspection cost as a percentage of production cost
	Mark-up for selling price (	) -must be greater than 1

following notations have been used in this study. Set-up time for a production cycle Idle time of a production cycle ̃ Fuzzy demand per year Holding cost per unit per year ($ per unit per year) Reliability of the production processwhich is known from the historical data of the production system Production lot size per normal cycle with reliability r Set up cost per cycle ($ per set-up) Production rate (units per year) in a 100% reliable system Production downtime for a cycle (set-up time + idle time) Number of cycles to recovery after the n th disruptiongiven from the management New disrupted cycle number from previous disruption Disruption period in the n th disruption Pre-disruption production quantity in the n th disruption Production time for Production quantity for a normal cycle i Production quantity for cycle i the recovery period after n th disruptionwhich is the decision variable; Productions up time for cycle a normal cycle i

Table 1 .

 1 1. Demand data as TFN and associated probabilities The problem is solved using the proposed GA based heuristic. The results are obtained from 30 different runs. The best recovery plan after single disruption is shown in Table2. The production system returns to original schedule from the sixth cycle after the disruption with , and so on. The maximum total profit in the recovery period is obtained as 1381112.5.

		Start
		Determine
		Assign		Randomly generate initial population of
				chromosomes
		Identify the first disrupted cycle
				Evaluate fitness function
		Initialize starting time from beginning
		of the disrupted cycle		Convergence
		Input cycle number with pre-disruption quantity and disruption period		GA operations: selection, crossover, mutation
		Solve the model to revise production		Update fitness function and individual
		quantity using GA
		Update	No	Meet stopping criteria?
				Yes
	Yes	Is there any other disruption?		Return the chromosome with revised produc-tion quantity
		No		Working principle of GA
		Stop

Table 2 .

 2 Best results obtained for the single disruption

	Disruption		Revised production quantity
	number						Total profit
	(n)
	1	5305	5638	6066	6469	6563	1381112.5

Table 3 .

 3 Data for the series of disruptions

	Disruption number (n)	Disrupted cycle number from previous disrup-tion	Pre-disruption quantity	Disruption peri-od
	1	1	1000	0.0045
	2	2	650	0.0092
	3	3	500	0.0025
	4	5	1500	0.0065
	5	2	0	0.0110
	6	4	800	0.0098
	7

Table 4 .

 4 Best results obtained for the series of disruptions

	Disruption		Revised production quantity
	number						Total profit
	(n)
	1	5585	6561	6545	6567	6563	1545915.7
	2	5272	5671	6111	6553	6572	1404616.5
	3	5611	6548	6570	6541	6550	1524257.8
	4	4804	6539	6433	6476	6522	1506722.6
	5	5271	5647	5994	6336	6381	1324912.7
	6	4941	5668	5936	6271	6532	1364154.2
	7	5657	5914	6049	6486	6443	1407321.9

Final Objective Function

Total profit, the objective function, is derived by subtracting all costs from the total revenues. Considering all the equations from ( 2) to [START_REF] Özekici | Inventory models with unreliable suppliers in a random environment[END_REF], the objective function is obtained as follows.

Max ̃ Total Revenues-Total Costs [START_REF] Mohebbi | An inventory model with non-resuming randomly interruptible lead time[END_REF] 4

Fuzzy Parameter

In this paper, we consider product demand as a triangular fuzzy number (TFN) to tackle uncertainty. A TFN ̃ is specified by a triplet ( ) and is defined by its continuous membership function ̃( ) ] as follows:

( ) and ( ) indicates the left and right branch of the TFN ̃ respectively. Ancut of ̃ can be expressed by the following interval [START_REF] Lee | Economic production quantity for fuzzy demand quantity, and fuzzy production quantity[END_REF]:

] The graded mean integration value (GMIV) of a LR-fuzzy number is introduced by Chen and Hsieh [START_REF] Chen | Graded mean integration representation of generalized fuzzy number[END_REF]. The graded mean integration representation method is based on the integral value of the graded mean -level of the LR-fuzzy number for defuzzifing LR-fuzzy numbers. By considering ̃ is a LR-fuzzy number and according to Chen and Hsieh [START_REF] Chen | Graded mean integration representation of generalized fuzzy number[END_REF], the GMIV of ̃ is defined as:

Disruption Recovery Model with Fuzzy Demand

In this section, fuzziness of demand is incorporated to the final mathematical model. The GMIV of the expected total profit function is evaluated. Relevant constraints are also developed with the GMIV of expected fuzzy demand. After simplifying the equation [START_REF] Mohebbi | An inventory model with non-resuming randomly interruptible lead time[END_REF], the following equation of the total profit is obtained:

Now, considering the fuzzy random demand ̃with the given set of ta ( ̃ ̃ ) ( ̃ ̃ ) ( ̃ ̃ ) ( ̃ ̃ ), the profit ( ̃) is also a fuzzy random variable and its expectation is a unique fuzzy number [START_REF] Bag | A production inventory model with fuzzy random demand and with flexibility and reliability considerations[END_REF] which is,

After integrating and simplifying the above equation of ( ̃), the GMIV of the total profit function, which is to be maximized, and obtained as:

Where,

GMIV of the expected fuzzy demand, ( ̃) Subject to the following constraints:

Solution Approach

We propose a genetic algorithm based heuristic to solve the model. Genetic algorithm is very popular technique to solve complex non-linear constrained optimization problem. GAs are general purpose optimization algorithms which apply the rules of natural genetics to explore a given search space [START_REF] Homaifar | Constrained optimization via genetic algorithms[END_REF]. The heuristic is designed to make a recovery plan from a single or a series of production disruptions. The proposed heuristic revises the production lot size of each cycle as long as disruptions take place in the system. For a series of disruptions, the heuristic revises the lot size of each cycle by considering the effect of all previous dependent disruptions. The proposed genetic algorithm based heuristic is presented in the Figure 2. The above mentioned heuristic is coded in MATLAB R2012a with the help of its optimization toolbox. In the proposed heuristic, following GA parameters are used to solve the model.