Ismael Bouassida Rodriguez

Karim Guennoun
email: guennoun@ehtp-ac.ma

Khalil Drira
email: khalil@laas.fr

Christophe Chassot
email: chassot@laas.fr

Mohamed Jmaiel
email: mohamed.jmaiel@enis.rnu.tn

Implementing a rule-driven approach for architectural self configuration in collaborative activities using a graph rewriting formalism

Keywords: D.2.11 [Software Engineering]: Software Architectures-Languages, D.2.9 [Software Engineering]: Management-Software configuration management, D.2.8 [Software Engineering]: Metrics-performance measures Autonomic computing, Rule-driven, Graph rewriting

Architectural adaptation is important for handling self-configuring properties of autonomic distributed systems. It can be achieved by model-based management of dynamic architectures. Describing dynamic architectures includes defining rules for reconfiguration management. Within this research context, several works have been conducted using formal specification to handle this complexity. Graph and graph rewriting-based approaches showed, through many studies, their appropriateness to tackle architectural adaptation problems. However, scalability of such approaches remains an open issue and has been rarely explored. In this paper, we investigate this issue. We introduce a graph-based general approach for handling of dynamic architectures, and we illustrate it within a scenario of collaboration support in Crisis Management Systems. We elaborate the formal models for dynamic architecture management. Using the French Grid GRID5000, we conducted an experimental study to assess the scalability of the elaborated models.

INTRODUCTION

The design of Self-configuring Systems that adapt their service to dynamically changing environments is among the main research directions of autonomic computing and communication [START_REF] Ganek | The dawning of the autonomic computing era[END_REF]. Providing such solutions for distributed software systems supporting group communication requires dynamically managing evolving group membership and dynamically connecting deployment nodes. It also requires dynamically distributing software services and components on such remotely interconnected deployment nodes.

For a number of group communication-based applications, to anticipate reconfiguration is important. In order to be applicable in different situations, one has also to ensure tractability of the elaborated solutions when changing from some to thousands users, nodes, components and services.

Providing generic and scalable solutions for automated self-reconfiguration in group collaboration support systems can be driven by rule-based reconfiguration policies. We elaborate a graph-based modelling approach and structural models that may represent the different interaction dependencies from different configuration-related point of views: communication flows between the distributed machines, the networked deployment nodes, and the service composition.

Our solution is based on graph grammar rewriting. Architectural reconfiguration is handled by defining architectures as graphs where vertices correspond to deployment nodes, software services and their internal components. It provides graph transformation to specify rules for changing deployment architecture while being in conformance to an architectural style. Dealing with dynamically evolving architectures requires at least describing the set of consistent architecture instances. This is mandatory to validate the management models and verify architectural constraint preservation. An architecture instance is considered consistent, if its corresponding description graph can be generated by a sequence of graph grammar productions. Our approach supports formal verification for correctness and safety proofs. We implemented a graph rewriting system that ensures automating our approach with a high performance making it tractable for large scale configurations.

To illustrate the proposed models and their transformations, we consider a case study of crisis management systems (CMS) involving several cooperating participants having different roles and functions. Related graph-based models are provided in section 4. Section 5 gives experimental results and analysis comments for scalability and performance of implemented models. Conclusion is provided in the last section.

RELATED WORK

Adaptation objectives, actions and properties are among the main facets of adaptability. They are studied and classified in this section.

Adaptability and adaptation

Adaptation is the operation to make changes to a program or an information system to maintain its functionalities and, if possible, to improve its performance in a certain execution environment. In the area of ubiquitous computing and context-awareness, adaptation is extended by the concept of adaptability that characterizes a system capability to change its behaviour to improve its performance or to continue its role in different environments. On the usefulness of adaptation, M. Satyanarayanan said in [START_REF] Satyanarayanan | Pervasive computing: Vision and challenges[END_REF] "adaptability is necessary when there is a significant disparity between supply and demand of a resource."

Goals

After developing an information system, many reasons can lead to adapt it. The reasons for this adaptation can be for corrective, evolutional or perfective purposes [START_REF] Ketfi | Adaptation dynamique, concepts et experimentations[END_REF]. In some cases, we can notice that the application does not behave properly or as expected. The corrective adaptation is a solution to identify the application module that causes the problem and to replace it by a new correct module. This new module provides the same functionality as the former. For the evolutional adaptation, when developing an application, some features are not taken into account. With the changing needs of the user, the application must be extended with new features. This extension can be accomplished by adding one or more modules to provide new features or modifying existing modules to enrich their functionality while keeping the same application architecture. The objective of the perfective adaptation is to improve application performance. For example, we can have a module that receives a lot of requests and fails to meet them. To avoid system performance degradation, we can duplicate this module to share the requests with the existing one.

Classification

The adaptation solutions suggested in the literature can be classified as follow.

Design vs run-Time

Two different adaptability views may be distinguished: the design time adaptability [START_REF] Dimka Karastoyanova | Extending web service flow models to provide for adaptability[END_REF][START_REF] Fahmy | Using graph rewriting to specify software architectural transformations[END_REF][START_REF] Ermel | Visual design of software architecture and evolution based on graph transformation[END_REF] and the run time adaptability [START_REF] Chang | Automatic configuration and run-time adaptation of distributed applications[END_REF][START_REF] Bade | Carsa -an architecture for the development of context adaptive retrieval systems[END_REF].

For the first view, we can find design support tools handles the application development cycle and optimizes the resources for example. For the run time adaptability [START_REF] Friday | Developing adaptive applications: The most experience[END_REF] presents several adaptation techniques which use proxy services, change model of interaction and reorganize application structure.

Local vs Distributed

Adaptation may have a local or a distributed scope. Adaptive components can be deployed on a single machine or distributed on several machines. In the first case, the adaptation is local and only local changes are performed. In the second case, it is distributed and synchronization problems between peer adaptive entities have to be managed [START_REF] Bridges | Supporting coordinated adaption in networked systems[END_REF].

Behavioural vs architectural

The adaptation solutions suggested in the literature distinguish behavioural and architectural aspects. The adaptation is behavioural (or algorithmic) when the behaviour of the adaptive service can be modified, without modifying its structure. Standard protocols such as TCP and specific protocols such as [START_REF] Wu | Streaming video over the internet: approaches and directions[END_REF][START_REF] Özgür | Atl: an adaptive transport layer suite for next-generation wireless internet[END_REF] provide behaviour-based adaptation mechanisms. Behavioural adaptation is easy to implement but limits the adaptability properties.

The adaptation is architectural when the service composition can be modified [START_REF]IEEE Std 1471-2000, IEEE Recommended practice for architectural description of software-intensive systems[END_REF][START_REF] Ellis | Toward a recommended practice for architectural description[END_REF] dynamically. In self-adaptive applications components are created and connected, or removed and disconnected during the execution. The architectural changes respond to constraints related to the execution context involving, for example, variations of communication networks and processing resources. They may also respond to requirement evolution in the supported activities involving, for example, mobility of users and cooperation structure modification.

Objects of adaptation

Adaptation approaches target many levels of information systems: User interfaces, Content, Services, Middleware and Transport.

User interface adaptation

User interface adaptation consist in producing Human-Computer means that can be deployed and used on different types of terminals and which meet the user's preferences. Major existing work in user interface adaptation is based on models that describe the different aspects of the interaction between humans and machines. These models are implemented in different XML or UML languages like UMLi [START_REF] Da Silva | UMLi: The unified modeling language for interactive applications[END_REF] and XIML [START_REF] Puerta | Ximl: a common representation for interaction data[END_REF]. These models are used to produce the adequate user interface code corresponding to the given XML or UML description. In the existing user interface adaptations, we distinguish two techniques: User Interface transformations and User Interface generation. In the first technique, the adaptation process starts from a description language which is very close to the user interface code that must be generated. This solution was adopted by [START_REF] Nasser | Adaptive bandwidth framework for provisioning connection-level qos for next-generation wireless cellular networks[END_REF] to produce adapted Web pages to different terminals starting from an XML description. In this kind of adaptation, style sheets are used to specify replacement rules of XML tags by scripts that can be directly used by the target device. The second approach consist in generating user interfaces code starting from a high level description which is completely independent from the target programming language of the user interface.

Content adaptation

Since the appearance of pervasive systems, the adaptation of multimedia content has been the subject of considerable research. Several techniques for adapting the delivered data to the user have been proposed. These techniques are based on textual transformations [START_REF] Nagao | Semantic annotation and transcoding: Making web content more accessible[END_REF], image transcoding [START_REF] Wee | Secure scalable streaming and secure transcoding with jpeg-2000[END_REF], or processing video and audio. One of the major issues in content adaptation is where the decision-making and transformations are made. In the literature, three general approaches have been proposed according to the location of adaptation processing between the source that hosts the content and destination that requests it: (i) on the content provider side, (ii) on the requester side, and (iii) at an intermediary (proxy) between the data source and the client.

The content provider-side solutions have some drawbacks. Indeed, the changes made on the content induce a calculation load and consequent resource consumption on the server. However, this approach is very suitable for situations with low variability and low adaptation frequency regarding the simplicity of its implementation. But it is not reliable for cases where the adaptation is triggered frequently.

The content requester-side approach is suitable when the transmission characteristics are less critical than the display limits of the user device [START_REF] Nandagopal | Service differentiation through end-to-end rate control in low bandwidth wireless packet networks[END_REF]. However, the usual complexity of adaptation processing hampers the wide adoption of this approach [START_REF] Perkis | Universal multimedia access from wired and wireless systems[END_REF]. Indeed, the client's terminal has generally very limited computing capacity, power and storage.

For the proxy solution, the flexibility of positioning the adaptation mechanisms on the best content distribution point is a major advantage compared to the other approaches (provider and content sides). However, the proxy must be a trusted party by the provider and the requester of the content. In addition, the third party may charge for the service it provides and the resources it employs to perform the adaptation for the receiver. Therefore, accounting mechanisms should be incorporated in the proxy solution in order to keep track of the amount of resources utilized and the usage of the data.

Service adaptation

Service-Oriented Architecture (SOA) paradigm is based on dynamically publishing and discovering services. This kind of architectures provides the possibilities to dynamically compose services for adapting applications to contexts. Service descriptions are published, via the registry, by service providers and dynamically discovered by service requesters. There are various implementation technologies like Web Services or COM/ DCOM (Component Object Model) [START_REF] Grimes | Professional Dcom Programming[END_REF] of Microsoft, the EJB (Enterprise Java Beans) [START_REF] Matena | Applying Enterprise Javabeans: Component-Based Development for the J2ee Platform[END_REF] of Sun Microsystems or CCM (CORBA Component Model) [START_REF] Omg | Corba components: Joint revised submission[END_REF] of OMG (Object Management Group).

Middleware adaptation

Other frameworks are proposed to provide adaptability for the middleware level. In [START_REF] Nasser | Adaptive bandwidth framework for provisioning connection-level qos for next-generation wireless cellular networks[END_REF], an adaptive framework for supporting multiple classes of multimedia services with different QoS requirements in wireless cellular networks is proposed. [START_REF] Sun | Cme: a middleware architecture for network-aware adaptive applications[END_REF] proposes CME, a middleware architecture for service adaptation based on network awareness. CME is structured as the software platform both to provide network awareness to applications and to manage network resources in an adaptive fashion.

Transport adaptation

At the transport level, [START_REF] Exposito | FPTP: the XQoS aware and fully programmable transport protocol[END_REF] provide frameworks for designing Transport protocols whose internal structure can be modified according to the application requirements and network constraints. Adaptation actions correspond to the replacement of a processing module or micro-protocol by another following a plug and play approach.

Model based adaptation approaches

There are many relevant contributions concerning system architecture adaptation. This kind of approaches uses model-based strategies to apply the necessary transformations on the systems architecture to adapt it to environment and requirement changes. These strategies define or reuse models describing the system software architectures. These models are also known as ADLs (Architecture description languages). We distinguish between there general ADL types: formal ADLs like graph grammars [START_REF] Hirsch | Graph grammars and constraint solving for software architecture styles[END_REF] and Petri nets [START_REF] Murata | Petri nets: Properties, analysis and applications[END_REF], semantic ADLs using ontologies [START_REF] Zhou | Applying ontology in architecture-based self-management applications[END_REF] and technical ADLs using XML deployment languages in general [START_REF] Dashofy | An infrastructure for the rapid development of xml-based architecture description languages[END_REF]. The technical ADLs can be proprietary or implementing the formal and the semantic architecture description models. These ADLs are used to guarantee the architecture evolving and correctness during the different predictable and unpredictable changes in the systems environment. The necessary actions to achieve such adaptations are specified in rules according to the application runtime context. [START_REF] Chaari | Adaptation in Context-Aware Pervasive Information Systems: The SECAS Project[END_REF] is an example of these approaches defining a complete model based architecture adaptation at the Service, content and user interface levels. [START_REF] Chassot | Towards autonomous management of qos through model-driven adaptability in communication-centric systems[END_REF] presents another model based method using graph grammars to adapt cooperative information systems to situation changes in the communication level.

Mobile networked systems for the support of group-wide activities have to provide dynamic adaptability for run-time and design-time changes. Such changes are induced by evolving requirements of the supported activities and by evolving resource constraints.

Designing and implementing self adaptive communicating systems to support such emerging activities is complex. Mastering this complexity may be achieved by adopting model-based design approaches associated with automated management techniques for dynamic architectural and behavioral adaptability management.

The works exposed allow clarifying the different problems related to self-adaptability. The proposed solutions have different features related to the targeted goals (e.g. QoS, security...), the considered layers (e.g. Application layer, Transport layer...), the adaptation actions (e.g. rate control, congestion control...), the adaptation properties (e.g. architectural, behavioural...), and the way to manage the adaptation and its autonomy.

To face the different kinds of requirements and constraints evolving, behavioral and architectural adaptability is required at several levels simultaneously. This typical cross layering problematic requires managing coordination without which it can lead to performances way below the targeted ones.

For group communication-based cooperative activities, ab-straction levels have to be identified. Adaptation at the highest levels should be guided by the evolving of the activity requirements. Adaptation at the lowest levels should be driven by the changes due to device/network constraints. Designing and implementing self-adaptive communicating systems is a complex task. To handle this complexity, several studies showed the need to lay on model-based design approaches associated with automated management techniques.

Static architectures are described by instances of components and interconnection links. The dynamic character of architectures requires additional description rules. Several works have addressed the dynamic architecture description, using different approaches [START_REF] Magee | Specifying Distributed Software Architectures[END_REF][START_REF] Allen | A formal basis for architectural connection[END_REF]. In order to guarantee the architecture evolving, correctness formal techniques are used. In particular, graphs represent a powerful expressive mean to specify respectively static and dynamic architecture aspects [START_REF] Metayer | Describing software architecture styles using graph grammars[END_REF][START_REF] Hirsch | Modeling software architectures and styles with graph grammars and constraint solving[END_REF]. For such approaches, graph vertices represent the software components, and the edges represent the links between these components. Dynamic architectures are described as graph grammars and architecture transformation is specified and ruled using graph rewriting models.

CASE STUDY

To expose the targeted problems and concepts, we consider the example of crisis management systems (CMS). In the sequel, we introduce this example and give two different execution steps and some related scenarios.

For CMS-like activities, collaboration is based on information exchange between mobile participants collaborating to achieve a common mission. A CMS team is composed of a number of participants having different roles: The controller of the mission, coordinators, and field investigators. Each group of investigators is supervised by a coordinator (Figure 1). Each participant is associated with an identifier, a role The controller's function consists in supervising the whole mission, i.e. deciding actions to be performed from the analysis of the observation feedbacks D transmitted by the coordinators. The scenario is divided into two steps: Investigation step and Action step. Initially, all investigator groups are in the "exploration step" where Investigators provide continuous feedbacks D to the coordinator; they also provide periodical feedbacks P and the coordinator sends continuous feedbacks P to the controller.

When an investigator meets a critical situation, the architecture of its related group has to be reconfigured to face this new situation. The group has to move to another execution step called an "action step" where The investigator that discovers the critical situation keeps sending both feedbacks D and P to the coordinator. It also provides feedbacks P to the other investigators of its group. Other investigators report feedbacks P to coordinator on the basis of feedbacks D transmitted by the critical investigator and the coordinator continue sending feedbacks P to the controller.

When the critical situation is resolved, the investigation group comes back to the exploration step.

GRAPH-BASED MODELS FOR ARCHI-TECTURAL ADAPTATION

In this section, we introduce our general framework of models. We use models based on graph grammar and graph rewriting to handle the reconfiguration management.

Graph grammars constitute a powerful and very expressive formalism for style description. Moreover, theoretical work on this field provides formal means to specify and check constraints on these architectures [START_REF] Rozenberg | Handbook of Graph Grammars and Computing by Graph Transformations[END_REF].Inspired from Chomsky's generative grammars [START_REF] Chomsky | Three models for the description of language. Information Theory[END_REF], graph grammars are defined, in general, as a classical system < AX; NT ; T ; P >, where AX is the axiom, NT is the set of the non-terminal vertices, T the set of terminal vertices, and P the set of transformation rules, also called grammar productions. An instance belonging to the graph grammar is a graph containing only terminal vertices and is obtained starting from axiom AX by applying a sequence of productions in P . There are different approaches for the definition of graph grammar production structure and the mechanisms used for the specification of their execution. In our approach,a graph grammar is of the form < AX; NT ; T ; P > where AX, N T and T keep the same significance as before. We use productions of type (L; K; R; C) where (L; K; R) corresponds to the structure of a Double PushOut (DPO) production [START_REF] Ehrig | Tutorial introduction to the algebraic approach of graph grammars[END_REF] and where C is a set of connection instructions. The instructions belonging to C are of the edNCE type [START_REF] Rozenberg | Handbook of Graph Grammars and Computing by Graph Transformations[END_REF]. They are specified by a system (n,δ,d,d ) where n corresponds to a vertex belonging to the daughter graph R, p and q are two edge labels, δ is a vertex label, and d and d are elements of the set in, out. For example, a production defined by the system (L; K; R; (n,δ,d,d )) is applicable to a graph G if it contains an occurrence of the mother graph L. The application of this production involves transforming G by deleting the subgraph (Del = L \ K) and adding the subgraph (Add = R \ K) while the subgraph K remains unchanged. All dangling edges will be removed.

The execution of the connection instruction implies the introduction of an edge between the vertex n belonging to the daughter graph R and all vertices n that are p-neighbours1 of and d-neighbours2 . This edge is introduced following the direction indicated by d and labelled by q.

An example is given in Figure 2. Production p3 has the same structure and transformation logic as in the DPO approach: node 5 is not removed even if it is matched with node 2' belonging to the L pattern because it also belongs to the K pattern. The example also considers the edNCE connection instructions c1, c2 and c3 allowing the correct addition of edges connecting nodes 2 and 6 and connecting nodes 3 and 7.

Following the commonly used conventions, we consider that vertices represent communicating entities (e.g. services, components) and edges correspond to their related interdependencies (e.g. communication links, composition dependencies).

An example of an architecture instance is given in figure 3. This instance includes, in addition to the controller (Cont), two coordinators (coord) that manage respectively three and two investigators (Inv). The graph edges are labelled by the exchanged data types (D/P ). Each participant have two attributes: the identifier and the deployment host. For example an investigator node is represented by a label like that inv(identif ier, deploymenthost). Architectural reconfiguration management consisting in transforming the system architecture evolving situation is formally ruled using graph rewriting techniques.

THE GRAPH GRAMMAR FOR ARCHI-TECTURE RECONFIGURATION

In the following, we provide an example of graph grammar for our case study to implement architecture reconfiguration.

Managing changes in activity organisation

The presented rewriting rule allows transforming the architecture from for moving the exploration step to the action step. The graph grammar is reduced to a single production grammar Pexp→act (table1) which is parameterized by the investigator identifier (here, noted A) who has discovered the critical situation. The architecture is transformed by splitting the communication channels between the coordinator and the other investigators into a communication channel of type P ,between these investigators and A into another communication channel of type D. Figure 4 gives an example of the application of this rule to move the architecture from the exploration step to the action step. Investigator i2 plays the role of the critical investigator A. Pexp→act allows the correct generation of the communication channels. The connection instruction ic1 allows to change the labels of the coordinator to the investigators from < D, P > to < P >. The connection instruction ic2 allows to connect the investigator A to the others investigator by an edge labelled (< P >).

Managing changes in group membership

We elaborate a graph grammar for our case study to manage the group extension. In the exploration step, the investigators produce D-data continuously and P-data periodically. The coordinators send P-data periodically to the controller. In the graph grammar GG (table2) we have four productions. ---→ Cont(c, m1), Inv(i, m3) The productions p1, p2 and p3 allow us to add new participants to the mission. The production p1 allows to initiate the group by connection an investigator, a coordinator and a controller. The production p2 allows to connect a coordinator and the edges necessary to that. According to the description of the case study, in the group, it is not possible to have a coordinator without an investigator. So the production p2 adds an investigator to the added coordinator. The production p3 allows to add an investigator to an existent coordinator. Such a propriety constitute the style of the architecture. The production p4 allows to finish the generation process and eliminate the non-terminal T emp. The example shown in figure 3, can be generated 3 by this sequence: p1;p2;p3;p3;p3;p4. The existence of this generation path prove the consistency of this architecture instance.

<D,P > -----→ Coord(co, m2)}) C = { }) p2 = (L = {Cont(c, m1), T emp}; K = {Cont(c, m1), T emp}; R = {Coord(co, m2), Inv(i, m3), Coord(co, m2) <P > ---→ Cont(c, m1), Inv(i, m3) <D,P > -----→ Coord(co, m2)}) C = { }) p3 = (L = {Coord(co, m1), T emp)}; K = {Coord(co, m1), T emp)}; R = {Inv(i, m2), Inv(i, m2) <D,P > -----→ Coord(co, m2) C = { }) p4 = (L = {T emp}; K = { }; R = { }) C = { })

EXPERIMENTAL STUDY

We conducted an experimental study of our graph rewriting system using the rules of the architecture adaptation models presented above.

These experiments have been achieved under grid GRID5000 [START_REF] Cappello | Grid'5000: a large scale, reconfigurable, controlable and monitorable Grid platform[END_REF], an experimental grid platform that interconnects several clusters geographically distributed in France. Considered grid nodes are single core machines (Sun Fire V20z) with an AMD Opteron 248 processor with 2.2 GHz of CPU speed 3 This is not the unique generation path and 2 GB of RAM. The provided network bandwidth is up to 10GB/s. We are studying scalability and limits of this model to investigate the tractability of an anticipation approach for adaptability.

The results, presented here, were obtained for graph grammar GG of section 4 and the reconfiguration rules executed to handle critical situation discovery (grammar production Pexp→act in table 1 of section 4). The experimental performance evaluation focuses execution time.

Experimentations for adapting architecture to organization changes

In this section we present, the results related to the execution of the Pexp→act rewriting rule. We are studying scalability and limits of our model-driven approach implementation.

We study the evolution of the execution time (figure 5), when the architecture description graph size growth. This corresponds to the application of production p1 of Pexp→act (Table 1). The implementation performs well and execution time is under 0.3 second for a 10 000 vertex graph (i.e. a system composed of 10 000 participants).

Experimentation adapting architecture to group membership changes

We generate the first N members of this graph grammar corresponding to the N first consistent instances of the ar- 2) whenever it is possible to obtain all consistent architecture graphs. The application of produc-p4 ensures the elimination of the non-terminal nodes in the graphs. The production applications include finding all graph matchings between the rules graph and transforming the architecture graph by adding nodes and edges. The production p2 adds two nodes and two edges. The production p1 adds one node and one edge. The system, takes about nine hours to calculate the 400 000 first consistent instances and more than 1 500 000 consistent instances are generated in about 194 hours. These values remains acceptable for design-time purposes. Using a unique grid node, the system generates almost the 400 000 first consistent instances in about 140 hours.

CONCLUSION

In this paper, we have studied a graph-based approach for implementing a model-based management of architectural reconfiguration. Our approach is useful when designing self-configuring software for large scale activities support. It supports scalable management of automated reconfiguration while allowing formal proofs to be developed for ensuring the consistency of the generated configurations. The models presented in this paper have been elaborated for the design of self-configuring Crisis Management Systems. The scalability study provided in this paper has been conducted to assess the tractability of our approach when dealing with large scale group applications.

Figure 1 :

 1 Figure 1: CMS architecture

Figure 2 :

 2 Figure 2: Combining DPO and edNCE

Figure 3 :

 3 Figure 3: Architecture instance (Exploration step)

Figure 4 :

 4 Figure 4: Reconfiguration from the exploration step to the action step

Table 1 :

 1 Pexp→act(A) = (p1, C = ic1, ic2) with:p1 = (L = {Coord(c, m1), Inv(A, m2), Inv(A, m2) <D,P > -----→ Cont(c, m1)}, K = { }, R = {Coord(c, m1), Inv(A, m2), Inv(A, m2) <D,P > -----→ Cont(c,m1)}) with:ic1 = (Coord(c, m1), < D, P > / < P >, Inv, in/in) ic2 = (Inv(A, m2), < D, P > / < D >, Inv, in/out) Production grammar Pexp→act GG = (AX, NT, T, P ) with: T = {Inv(i, m3), Coord(co, m), Cont(c, m)}, and NT = {T emp}, and P = {p1, . . . , p4} p1 = (L = {AX}; K = { }; R = {Cont(c, m1), Coord(co, m2), Inv(A, m3), T emp), Coord(co, m2) <P >

Figure 5 :

 5 Figure 5: Execution time of applying p1 of Pexp→act wrt architecture graph size

Figure 6 :

 6 Figure 6: Execution time of applying p1, p2, p3 and p4 of GG wrt graph number

Table 2 : Group membership change grammar GG

 2

p-neighbours of a vertex n are all vertices n such that there exists an edge labelled by p which connects n and n .

In-neighbours if d = in and out-neighbours otherwise.


Souad Bezzaoucha

Benoît Marx

Didier Maquin

José Ragot

Stabilization of Nonlinear Systems Subject to Actuator Saturation

Keywords: Takagi-Sugeno systems, actuator saturation, parallel distributed compensation control law

This paper addresses the stabilization of nonlinear systems described by Takagi-Sugeno models affected by input actuator saturation. A parallel distributed compensation design is used for the state feedback controller. Stabilization conditions in the sense of the Lyapunov method are derived and expressed as a linear matrix inequality problem. The obtained gains depend on the actuator saturation limits. A cart-pendulum example is presented to illustrate the effectiveness of the proposed approach.

I. INTRODUCTION

Actuator saturation or control input saturation is probably one of the most usual nonlinearity encountered in control engineering due to the physical impossibility of applying unlimited control signals and/or safety constraints. In general, there are two main design strategies to deal with actuator saturations. The first strategy is a two-step approach in which a nominal linear controller is first constructed by ignoring the actuator saturation. Then, after this controller has been designed, usually using standard linear design tools, a so called anti-windup compensator is designed to handle the saturation constraints [START_REF] Syaichu-Rohman | Anti-windup schemes for discrete time systems: An LMI-based design[END_REF], [START_REF] Zaccarian | A common framework for anti-windup, bumpless transfer and reliable design[END_REF], [START_REF] Mulder | Multivariable antiwindup controller synthesis using linear matrix inequalities[END_REF], [START_REF] Grimm | Antiwindup for stable linear systems with input saturation: an LMI-based synthesis[END_REF]. A typical anti-windup scheme consists in augmenting a nominal pre-designed linear controller with a compensator based on the discrepancy between unsaturated and saturated control signals fed to the plant [START_REF] Syaichu-Rohman | Optimisation based feedback control of input constrained linear systems[END_REF]. The second strategy considers the saturation constraint from the beginning of the design task. Several approaches were developed, one of them is the invariant sets framework, which has been significantly developed in control engineering over the last decades, see [START_REF] Alessandri | On estimation error bounds for receding-horizon filters using quadratic boundedness[END_REF], [START_REF] Blanchini | Set-theoretic methods in control. Systems and Control: Foundations and Applications[END_REF]. This framework ensures that every state trajectory that starts inside the invariant set will not exceed it, i.e. the system states will be bounded inside this set. An interesting approach for the determination of the invariant sets framework was applied in [START_REF] Cao | Robust stability analysis and fuzzy-scheduling control for nonlinear systems subject to actuator saturation[END_REF] and [START_REF] Tarbouriech | Stability and Stabilization of Linear Systems with Saturating Actuators[END_REF]. This approach consists in the so-called polytopic rewriting of the saturation constraint. This polytopic representation is then used to determinate (maximize) the largest invariant set in which the system states remain bounded. In the contrary of the previous method, that most likely ensures to never reach the saturation level in order to preserve the closed loop performances, in this paper, we explicitly consider the saturation constraint, and even admit its occurrence. In fact, the authors use the Takagi-Sugeno (T-S) representation of the saturation (as well known as the polytopic representation) to integrate the limitation constraints into the control synthesis, such that the system stability is ensured and the control gains are calculated depending on the saturation level. The authors would like to precise that even if the expressions of the T-S saturation has some similarity to the polytopic one used in [START_REF] Cao | Robust stability analysis and fuzzy-scheduling control for nonlinear systems subject to actuator saturation[END_REF] and [START_REF] Tarbouriech | Stability and Stabilization of Linear Systems with Saturating Actuators[END_REF], the development, control strategy and objectives are completely different since in the proposed approach the invariant sets are never considered; moreover the objective is the stabilization of nonlinear systems represented with the Takagi-Sugeno (T-S) models and the synthesis of a state feedback controller by parallel distributed compensation (PDC) with control gains explicitly depending on the saturation level. The controller is subject to actuator saturation and the input saturation is directly taken into account in the controller design. Stabilization conditions are derived from the Lyapunov method and expressed as linear matrix inequalities (LMI). A cart-pendulum example is considered to illustrate the effectiveness of the proposed approach. The nonlinear equations of the cart-pendulum are given and, using the Sector Nonlinearity Transformation (SNT) [START_REF] Tanaka | Fuzzy Control Systems Design and Analysis: A Linear Matrix Inequality Approach[END_REF], [START_REF] Nagy | Systematic multimodeling methodology applied to an activated sludge reactor model[END_REF], a T-S model of the system is deduced and used. The system is subject to actuator saturations, it will be shown that these input constraints may cause the instability of the nonlinear system, but with the proposed approach the stability of the closed-loop system is ensured. The rest of this paper is organized as follows. Section 2 introduces the Takagi-Sugeno structure for modeling and some preliminary results, mathematical notations and a brief description of the saturation. It is followed by the representation of the nonlinear saturation by a T-S structure in section 3. In section 4 is designed a state feedback control law depending on the saturation bounds. A numerical example and some simulation results are given in section 5. Conclusions and future works are exposed in section 6.

II. PRELIMINARIES

A. Takagi-Sugeno structure for modeling Initially introduced in [START_REF] Takagi | Fuzzy identification of systems and its applications to modeling and control[END_REF], the T-S models represent a simple and accurate method to study the nonlinear behaviors. The T-S modeling allows to represent the behavior of nonlinear systems by the interpolation of a set of linear submodels [START_REF] Tanaka | Fuzzy Control Systems Design and Analysis: A Linear Matrix Inequality Approach[END_REF], [START_REF] Tanaka | A multiple Lyapunov function approach to stabilization of fuzzy control systems[END_REF], [START_REF] Tanaka | A descriptor system approach to fuzzy control system design via fuzzy Lyapunov functions[END_REF]. Each submodel contributes to the global behavior of the nonlinear system through a weighting function µ i (ξ(t)).

The T-S structure is given by

         ẋ(t) = n i=1 µ i (ξ(t))(A i x(t) + B i u(t)) y(t) = n i=1 µ i (ξ(t))(C i x(t) + D i u(t)) (1)
where x(t) ∈ R nx is the system state variable, u(t) ∈ R nu is the control input and y(t) ∈ R m is the system output. ξ(t) ∈ R q is the decision variable vector assumed to be measurable (as the system output) or known (as the system input). The weighting functions µ i (ξ(t)) of the n submodels satisfy the convex sum property

     n i=1 µ i (ξ(t)) = 1 0 ≤ µ i (ξ(t)) ≤ 1, i = 1, . . . , n (2)
In the remaining of the paper, the following lemmas are used: Lemma 1. Consider two matrices X and Y with appropriate dimensions and G a symmetric positive definite matrix. The following property is verified

X T Y + Y T X ≤ X T GX + Y T G -1 Y (3)
Lemma 2. (Congruence) Consider two matrices X and Y , if X is positive (resp. negative) definite and if Y is a full column rank matrix, then the matrix Y XY T is positive (resp. negative) definite.

B. Mathematical notations

The following notations are used throughout the paper: a block diagonal matrix with the square matrices A 1 , . . . , A n on its diagonal is denoted diag(A 1 , . . . , A n ). For any matrix, M , S(M ) is defined by S(M ) = M + M T .The smallest and largest eigenvalues of the matrix M are respectively denoted λ min (M ) and λ max (M ). The saturation function for a signal ν(t) is defined by [START_REF] Grimm | Antiwindup for stable linear systems with input saturation: an LMI-based synthesis[END_REF], where ν max and ν min denote the saturation levels.

sat(ν(t)) :=    ν(t) if ν min ≤ ν(t) ≤ ν max ν max if ν(t) > ν max ν min if ν(t) < ν min (4)

III. PROBLEM STATEMENT

A. Takagi-Sugeno saturation control

The main idea of this work is to model the nonlinear actuator saturation using the Takagi-Sugeno representation (section II-A) and then propose a PDC control law ensuring stability of the closed loop system. For that, it is proposed to re-write the saturation equation ( 4) for each component of the control input vector under a particular form. Let us consider a control input vector u(t) ∈ R nu , defined by

u(t) = ( u 1 (t) . . . u nu (t) ) T (5)
The control input under actuator saturation constraint is given by

u sat (t) = u 1 sat (t) . . . u nu sat (t) T ( 6
)
where u j sat = sat(u j (t)), for j = 1, . . . , n u . Considering the three parts of the saturated signal (4), each component of the vector u sat (t) is written as

u j sat (t) = 3 i=1 µ j i (u j (t)) (λ j i u j (t) + γ j i ), j = 1, . . . , n u (7) with λ j 1 = 0 λ j 2 = 1 λ j 3 = 0 (8) γ j 1 = u j min γ j 2 = 0 γ j 3 = u j max (9)
and the weighting functions defined by

       µ 1 (u j (t)) = 1-sign(uj (t)-u j min ) 2 µ 2 (u j (t)) = sign(uj (t)-u j min )-sign(uj (t)-u j max ) 2 µ 3 (u j (t)) = 1+sign(uj (t)-u j max ) 2 (10)
Based on the convex sum property of the weighting functions (2), the control input vector u(t) ∈ R nu subject to actuator saturation can be written in order to have the same activation functions for all the input vector components as:

u sat (t) =                      3 i=1 µ 1 i (λ 1 i u 1 (t) + γ 1 i ) ×   nu k=2 3 j=1 µ k j (u k (t))   . . . 3 i=1 µ i (λ i u (t) + γ i ) ×   nu k=1,k = 3 j=1 µ k j u k (t)   . . . 3 i=1 µ nu i (λ nu i u nu (t) + γ nu i ) ×   nu-1 k=1 3 j=1 µ k j u k (t)                        (11)
For n u inputs, 3 nu submodels are obtained. Thus, it is important to note that ( 11) is an analytical expression of the actuators saturation directly expressed in term of the control variable, that can also be expressed as

u sat (t) = 3 nu i=1 µ i (t)(Λ i u(t) + Γ i ) (12)
The global weighting functions µ i (t), the matrices

Λ i ∈ R nu×nu and vectors Γ i ∈ R nu×1 are defined as follows              µ i (t) = nu j=1 µ j σ j i (u j (t)) Λ i = diag(λ j σ j i ) Γ i = γ 1 σ 1 i , . . . , γ nu σ nu i T (13)
where the indexes σ j i (i = 1, . . . , 3 nu and j = 1, . . . , n u ), equal to 1, 2 or 3, indicate which partition of the j th input (µ j 1 , µ j 2 or µ j 3 ) is involved in the i th submodel. The relations between i the number of the submodel) and the σ j i indices are given by the following equation

i = 3 nu-1 σ 1 i +3 nu-2 σ 2 i +. . .+3 0 σ nu i -(3 1 +3 2 +. . .+3 nu-1 )
The indices σ j i are such that ((σ 1 i -1), . . . , (σ nu i -1)) correspond to (i-1) in base 3. An illustrative example is given for two inputs (n u = 2), with

u sat (t) = u 1 sat (t) u 2 sat (t) T (14)
Since three partitions are defined for each input, the Takagi-Sugeno model for u sat (t) is then represented by 3 2 submodels

u sat (t) = 9 i=1 µ i (t)(Λ i u(t) + Γ i ) (15)
with the parameters µ i , Λ i and Γ i given by the following table

submodel i (σ 1 i , σ 2 i ) µi(t) Λi Γi 1 (1, 1) µ 1 1 µ 2 1 diag(λ 1 1 , λ 2 1 ) γ 1 1 γ 2 1 T 2 (1, 2) µ 1 1 µ 2 2 diag(λ 1 1 , λ 2 2 ) γ 1 1 γ 2 2 T 3 (1, 3) µ 1 1 µ 2 3 diag(λ 1 1 , λ 2 3 ) γ 1 1 γ 2 3 T 4 (2, 1) µ 1 2 µ 2 1 diag(λ 1 2 , λ 2 1 ) γ 1 2 γ 2 1 T 5 (2, 2) µ 1 2 µ 2 2 diag(λ 1 2 , λ 2 2 ) γ 1 2 γ 2 2 T 6 (2, 3) µ 1 2 µ 2 3 diag(λ 1 2 , λ 2 3 ) γ 1 2 γ 2 3 T 7 (3, 1) µ 1 3 µ 2 1 diag(λ 1 3 , λ 2 1 ) γ 1 3 γ 2 1 T 8 (3, 2) µ 1 3 µ 2 2 diag(λ 1 3 , λ 2 2 ) γ 1 3 γ 2 2 T 9 (3, 3) µ 1 3 µ 2 3 diag(λ 1 3 , λ 2 3 ) γ 1 3 γ 3 3 T TABLE I. WEIGHTING FUNCTIONS, MATRICES Λ i AND Γ i , FOR nu = 2

B. Problem statement

Let us now consider a T-S nonlinear system represented by the following state equation

ẋ(t) = n i=1 µ i (ξ(t))(A i x(t) + B i u(t)) (16)
The control input u(t) is subject to actuator saturation, then the system (16) becomes

ẋ(t) = n i=1 µ i (ξ(t))(A i x(t) + B i u sat (t)) (17)
From [START_REF] Takagi | Fuzzy identification of systems and its applications to modeling and control[END_REF], equation (17) can be written as

ẋ(t) = n i=1 3 nu k=1 µ i (ξ(t))µ sat k (t)(A i x(t) + B i (Λ k u(t) + Γ k ))
(18) with n the number of sub-models, n u the dimension of the input vectors, µ i (ξ) and ξ respectively are the system weighting functions and the measurable decision variables. µ sat k , Λ k and Γ k are respectively the saturation weighting functions and the parameters defined by equations [START_REF] Tanaka | A multiple Lyapunov function approach to stabilization of fuzzy control systems[END_REF].

IV. SATURATED STATE FEEDBACK CONTROL INPUT

The objective is to design a stabilizing time-varying state feedback controller ensuring the stability of the system, even in the presence of control input saturation. The solution is obtained by representing the saturation as a T-S system and by solving an optimization problem under LMI constraints. In the nominal case, if no saturation is affecting the system (1), the well-known following PDC (Parallel distributed compensation) can be applied

u(t) = - n j=1 µ j (ξ(t))K j x(t) (19)
Using the quadrative Lyapunov function

V (x(t)) = x T (t)P x(t), P = P T > 0 ( 20
)
the stability of the closed-loop system is obtained by computing the gains K j from [10]

P 1 A T i +A i P 1 -R T j B T i -B i R j < 0 i = 1, . . . , n; j = 1, . . . , n (21)
with K j = P -1 1 R j (where P 1 = P -1 ). In the presence of input saturation, the gains computed from (21) do not ensure the closed-loop stability. The objective is then to design a nonlinear state feedback controller (19) in order to guarantee the stability of the saturated system (18) such that the control gains K i depend on the saturation limits. By replacing the control law (19) in the T-S system (18), the obtained closed-loop system is the following

ẋ(t) = n i=1 n j=1 3 nu k=1 µ i (ξ(t))µ j (ξ(t))µ sat k (t) ((A i -B i Λ k K j )x(t) + B i Γ k ) (22)
The determination of the controller gains is now detailed.

Theorem 1. There exists a time-varying state feedback controller (19) for a saturated input system (18) ensuring that the system state converges toward an origin-centred ball of radius bounded by β if there exists

P 1 = P T 1 > 0, R, Σ k = Σ T k > 0 solutions of the following optimization problem min P1, R, Σ k β (23) s.t. Q ijk I I -βI < 0 (24)
with

Q ijk = S(A i P 1 -B i Λ k R j ) I I -Σ k (25) and Γ T k B T i Σ k B i Γ k < β (26)
for i = 1, . . . , n, j = 1, . . . , n, k = 1, . . . , 3 nu . The controller gains are given by

K j = P -1 1 R j (27)
Proof:

Let us define the quadratic Lyapunov function (20). According to ( 22) and (20), the time derivative of V (x(t)) is given by

V (x(t)) = n i=1 n j=1 3 nu k=1 µ i (ξ(t))µ j (ξ(t))µ sat k (t) S(x T (t)P B i Γ k + x T (t)P (A i -B i Λ k K j )x(t)) (28)
Using Lemma 1, it follows that

S(x T (t)P B i Γ k ) ≤ Γ T k B T i Σ k B i Γ k + x T (t)P Σ -1 k P x(t) (29)
The time derivative of the Lyapunov function ( 28) is bounded as follows

V (x(t)) ≤ n i=1 n j=1 3 nu k=1 µ i (ξ(t))µ j (ξ(t))µ sat k (t) x T (t)(S(P (A i -B i Λ k K j )+P Σ -1 k P )x(t)+Γ T k B T i Σ k B i Γ k (30)
Let us define

Q ijk = S(P (A i -B i Λ k K j )) + P Σ -1 k P (31) ε = min i=1:n,j=1:n,k=1:3 nu λ min (-Q ijk ) (32) δ = max i=1:n,k=1:3 nu Γ T k B T i Σ k B i Γ k (33)
Since Σ k > 0 and from equation ( 30), according to Lyapunov stability theory [START_REF] Zhang | A new approach to observer-based fault-tolerant controller design for Takagi-Sugeno fuzzy systems with state delay[END_REF],

V (t) < -ε x 2 +δ. It follows that V (t) < 0 for    Q ijk < 0 and x 2 > δ ε (34)
which means that x(t) is uniformly bounded and converges to a small origin-centered ball of radius δ ε . Applying Lemma 2, Q ijk < 0 is equivalent to

S(P -1 A T i -P -1 K T j Λ T k B T i ) + Σ -1 k < 0 (35)
Using the variable changes

P 1 = P -1 R j = K j P 1 (36)
condition ( 35) is linearized as

S(P 1 A T i -R T j Λ T k B T i ) + Σ -1 k < 0 i = 1, . . . , n; j = 1, . . . , n; k = 1, . . . , 3 nu (37)
Applying Schur's complement to (37), it becomes

Q ijk = S(P 1 A T i -R T j Λ T k B T i ) I I -Σ k < 0 (38)
As the weighting functions satisfy (2) and Σ k > 0, if (38) is satisfied for i = 1, . . . , n, j = 1, . . . , n, k = 1, . . . , 3 nu , and x 2 > δ ε , then V (x(t)) < 0, which implies that x(t) converges to an origin centered ball of radius δ ε . The objective is now to minimize the radius δ ε . Firstly δ is bounded by β (33) and the LMIs (26). From (24), with a Schur complement, it obviously follows that

(1/β) I < -Q ijk , i = 1, . . . , 3 nu (39)
implying that all the eigenvalues of (-Q ijk ) are larger that 1/β. As a consequence 1/β < ε holds, and finally the radius is bounded by β.

Remark 1. It is important to highlight that the proposed approach ensures the stability of nonlinear systems that may be destabilized by saturated control with a control proportional to the saturation limits. If the submodels are initially unstable, the proposed approach is not suitable since the LMI condition cannot be fulfilled. Thus, to overcome the submodels stability constraint, a solution would be to consider a dynamic output feedback control instead of the state feedback one [START_REF] Bezzaoucha | Contribution to the constrained output feedback control[END_REF].

V. NUMERICAL EXAMPLE

Let us consider a nonlinear cart-pendulum system illustrated in figure 1. The pendulum rotates in a vertical plan around an axis located on a cart. The cart can move along a horizontal rail, lying in the rotation plane. The state of the system is a vector x = ( x 1 x 2 x 3 x 4 )

T where x 1 = z(t) is the cart position, x 2 = ż(t) the cart velocity, x 3 = θ(t) is the angle between the downward direction and the pendulum (x 2 = 0 for the upright position of the pendulum), and x 4 = θ(t) is the pendulum angular velocity. A control force F (t) parallel to the rail is applied to the cart. The pendulum mass and cart mass are denoted m = 1Kg and M = 5Kg respectively. l = 0.1m is the pendulum length and I the moment of inertia of the pendulum with respect to its axis on the cart. The cart friction is composed of the viscous friction proportional to the cart velocity f ż(t) (f = 100N/m/s), a friction proportional to the cart position k s z(t) (k s = 0.001N/rad/s) and a friction torque in the angular motion of the pendulum proportional to the angular velocity k θ(t) (k = 0.003N/rad/s). The state equations are as follows:

               (m + M )z(t) + k s z(t) + f ż(t) -ml θ(t)cos(θ(t)) +ml θ2 (t)sin(θ(t)) = F (t)
-mlz(t)cos(θ(t)) + (ml 2 + I) θ(t) +k θ(t) + mglsin(θ(t)) = 0 (40) In order to ensure the nonlinear system stability, the system equations ( 40) is written under a T-S from by applying the sector nonlinearity approach which allows to exactly represent the nonlinear system without any loss of information. Since the transition from the nonlinear system to the T-S one is not the main subject of this paper and for place limitation, only the main steps are given. First, the quasi-LPV form (quasi-linear with parameter varying) is derived

ẋ(t) = A(x(t))x(t) + B(x(t))u(t) (41)
with :

x(t) = z(t) ż(t) θ(t) θ(t) T , u(t) = F (t) (42) A(x(t)) =     0 1 0 0 -k s /(m + M ) a 1 (t) a 2 (t) a 3 (t) 0 0 0 1 a 7 (t) a 4 (t) a 5 (t) a 6 (t)     B(t) = ( 0 b 1 (t) 0 b 2 (t) )
T (43) and the parameters a i (t) and b i (t) are defined as follows:

                                       a 1 (t) = -f m+M - mlf z1(t)z 2 2 (t) m+M a 2 (t) = -mlz3(t)z4(t) m+M - m 2 l 2 z 2 (t)z1(t)z3(t)z4(t) m+M -mlgz 1 (t)z 2 (t)z 3 (t) a 3 (t) = kz 1 (t) a 4 (t) = -f z 1 (t)z 2 (t) a 5 (t) = -mlz 1 (t)z 2 (t)z 3 (t)z 4 (t) -(m + M )gz 1 (t)z 3 (t) a 6 (t) = -k(m+M ) ml z 1 (t) a 7 (t) = -k s z 1 (t)z 2 (t) b 1 (t) = 1 m+M + ml m+M z 2 2 (t)z 1 (t) b 2 (t) = z 1 (t)z2(t)
(44) with the following non-linearity :

         z 1 (t) = 1 -mlcos 2 (x3(t))+(l+I/ml)(m+M ) z 2 (t) = cos(x 3 (t)) z 3 (t) = sin(x3(t)) x3 z 4 (t) = x 2 4 (t) (45)
The reader will note that the nonlinearities choice (z 1 (t), z 2 (t), z 1 (t) and z 4 (t)) is not unique and that the premise variables are bounded and given by the angle and speed limitation. Now, the convex polytopic transformation is applied

z i (t) = F i1 (t)z i1 + F i2 (t)z i2 , i = 1, . . . , 4 (46)
with

z i1 = max z i (t), z i2 = min z i (t), i = 1, . . . , 4 (47)
F i1 (t) = zi(t)-zi2 zi1-zi2 F i2 (t) = zi1-zi(t) zi1-zi2 (48)
The two partitions F i1 (t) and F i2 (t) of each premise variable z i (t), i = 1, 2, 3, 4 will contribute to the reconstruction of each weighting function of the T-S model. The sub-model number depends on the non-linearities, such that r = 2 4 = 16.

The T-S model is then given by the following equations:

ẋ(t) = 16 i=1 µ i (t)(A i x(t) + B i u(t)) (49) with          µ i (t) = 16 j=1 F j,σ j i (z j (t)) A i = A(z 1,σ 1 i , z 2,σ 2 i , z 3,σ 3 i , z 4,σ 4 i ) B i = B(z 1,σ 1 i , z 2,σ 2 i , z 3,σ 3 i , z 4,σ 4 i ) (50)
σ j i may take the value 1 or 2 and represents the number of the j th partition. In order to simplify the LMI implementation, it is assumed that the control gains K j are equal for all the submodels. Then, two control laws are proposed:

u(t) = -K n x(t) and u(t) = -K T S x(t) (51)
For an initial position x 0 = ( 0 0 π/12 0 ) T , figure 2 shows the system states for the following three cases: x n for the nominal case without saturation, x nsat is the system states with saturated nominal control, it is clear from the depicted figure that the saturation has a destabilizing effect on the system. And finally x satT S denotes the system state obtained with the proposed T-S approach. From the depicted figures, it is clear that with the proposed T-S approach the system stability is ensured as well as the convergence to an origin centered-ball. However, since we considered the same gain for all the 16 sub-models, the obtained solution is slightly conservative, which may justify the oscillatory response for the third and fourth state. In order to relax these constraints, the nonlinear system may be simplified by considering the small angles approximations: cos(θ) ≈ 1 and sin(θ) ≈ θ. Hence, instead of the four previous non-linearities, only one is considered, namely: z 4 (t) = x 2 4 (t) and a T-S system with two sub-models is then defined with the T-S controller is then given by u(t) = -

As seen on Figure 3, displaying the obtained state responses, x 4nsat

x 4satTS Fig. 3. System states it is clear that the system states are stable with less oscillations and converge to an origin centered ball.

VI. CONCLUSIONS

Considering the saturation nonlinearity, a nonlinear system can be represented in an augmented T-S form. It is important to note that, using the proposed representation of the actuator saturation, these saturations are expressed in terms of the control variables and of the saturation limits. As a consequence, these limits are taken into account when solving the LMI optimization problem and when computing the gains of the controller. Moreover, the proposed Takagi-Sugeno approach allows to extend the use of linear tools, namely the LMI formalism, to a nonlinear control problem. As an illustrative example, the state feedback controller design was applied to a cart-pendulum nonlinear T-S system. The obtained results show the effectiveness of the proposed approach for the studied example.

Fig. 1 .

 1 Fig. 1. Cart-pendulum system

3 ,

 3 for the nominal controller (without input saturation) and the proposed T-S controller respectivelly. The nominal control gain calculated from (21) is equal to: K n = ( 3.4068 -95.7038 -29.8001 3.8235 ) (52) Considering the saturation levels equal u min = 0 and u max = The control gain K T S is then calculated by solving the LMI given by theorem 1 and is equal to: K T S = ( 0.0003 5.9835 0.6139 -0.0024 )

2 i=1µFig. 2 .

 22 Fig. 2. System states


domain: phys spi

text:
Clément Bielinski

Badr Kaoui
email: badr.kaoui@utc.fr

Numerical method to characterise capsule membrane permeability for controlled drug delivery

Keywords: Controlled release, absorption, capsules, membrane permeability, computer simulations

Design and characterisation of capsules is not an easy task owing to the multiple involved preparation factors and parameters. Here, a novel method to characterise capsule membrane permeability to solute molecules by an inverse approach is proposed. Transport of chemical species between the capsule core and the surrounding medium through the membrane is described by the Fick's second law with a position-dependent diffusion coefficient. Solutions are computed in spherical coordinates using a finite difference scheme developed for diffusion in multilayer configuration. They are validated using semi-analytical solutions and fully three-dimensional lattice Boltzmann simulations. As a proof of concept, the method is applied to experimental data available in the literature on the kinetics of glucose release and absorption to determine the membrane permeability of capsules. The proposed method is easy to use and determines correctly the permeability of capsule membranes for controlled drug release and absorption applications.

I. INTRODUCTION

The last decades have known an increasing interest in using capsules, closed polymeric membranes, to encapsulate active agents, which has led to a large amount of publications and patents. Capsules are composed of an inner core, where active molecules are encapsulated, and an outer protective thick shell or thin membrane, as illustrated in Fig. 1. The main advantage of capsules is their ability to monitor the release of their cargo in a controlled manner. Moreover, their membrane improves the mechanical robustness and protects the encapsulated material against undesired external chemical and mechanical damages. Applications of capsules span from food industry to medicine, with capsule size ranging from few millimeters down to few nanometers. They are fabricated through diverse techniques [START_REF] Ghosh | Functional Coatings[END_REF][START_REF] Sagis | Microencapsulation and Microspheres for Food Applications[END_REF], including emulsification in flow focusing microfluidic devices [START_REF] Chu | Fabrication and in situ characterization of microcapsules in a microfluidic system[END_REF][START_REF] Galle | On-chip encapsulation via chaotic mixing[END_REF][START_REF] Li | Microfluidic fabrication of microparticles for biomedical applications[END_REF]. Capsules with more complex structures are being developed by adding multiple extra layers or chemically reinforced shells to meet specific performance requirements [START_REF] Daehne | Fabrication of micro reaction cages with tailored properties[END_REF][START_REF] Qiu | Studies on the drug release properties of polysaccharide multilayers encapsulated ibuprofen microparticles[END_REF][START_REF] Zarket | Onion-like multilayered polymer capsules synthesized by a bioinspired inside-out technique[END_REF]. Capsules are also used as miniaturised bioreactors, where living cells are enclosed to carry biochemical reactions as a response to specific applied stimuli [START_REF] Ma | Core-shell hydrogel microcapsules for improved islets encapsulation[END_REF][START_REF] Fischer | Triggered release of loads from microcapsule-in-microcapsule hydrogel microcarriers: en-route to an "artificial pancreas[END_REF]. Scheme showing the structure of a core-shell capsule with radius R. Rc is the radius of the core and R b the radius of the bulk where the capsule is immersed. The shell is supposed to be homogeneous with a uniform thickness δ. The scheme also shows the mass transport directions of the release and absorption of a solute.

Design and optimisation of capsules is challenging due to the multiple physico-chemical and geometrical parameters involved. The release kinetics of the encapsulated solute is influenced by the capsule membrane properties, which have to be optimised in order to deliver molecules at desired rates. For example, for medical applications, a poorly designed capsule may lead to an inefficient treatment, and even worse, to drug concentration exceeding the toxicity threshold. The main role of capsule membrane is to moderate the exchange of chemical species between the capsule inner core and the outer surrounding medium. This depends strongly on the permeability of the membrane, which is paradoxically not known a priori. It is dictated by multiple factors, such as the local curvature, the porosity and the tortuosity of the membrane resulting from the quality of the cross-linking of the polymer, of which the membrane is made. Thus, the permeability is not a universal determined property. It depends rather on the material and the conditions under which a capsule is fabricated. It differs from a population of capsules to another due to slight accidental changes in the fabrication process. The capsule permeability cannot be measured directly because of the impossibility of accessing the protected inner core, for example, to place a probe apparatus. The closed topological character of the capsules leaves only the option to access physically their external surrounding environment.

The state-of-the-art approach to determine capsule permeability consists in measuring the evolution over time of the solute concentration in the external surrounding medium, which is then fit with an equation. This approach has been employed in many studies [START_REF] Dembczynski | Characterization of small molecules diffusion in hydrogel-membrane liquid-core capsules[END_REF][START_REF] Koyama | Evaluation of mass-transfer characteristics in alginate-membrane liquid-core capsules prepared using polyethylene glycol[END_REF][START_REF] Yao | Mass transfer behavior of solutes in NaCS-PDMDAAC capsules[END_REF][START_REF] Ameloot | Interfacial synthesis of hollow metal-organic framework capsules demonstrating selective permeability[END_REF][START_REF] Henning | Sealing liquid-filled pectinate capsules with a shellac coating[END_REF][START_REF] Leick | Alginate/poly-L-lysine capsules: mechanical properties and drug release characteristics[END_REF][START_REF] Shahravan | Encapsulation and controlled release from core-shell nanoparticles fabricated by plasma polymerization[END_REF][START_REF] Rolland | Physicochemical properties of aqueous core hydrogel capsules[END_REF][START_REF] Ben Messaoud | Physico-chemical properties of alginate/shellac aqueous-corecapsules: Influence of membrane architecture on riboflavin release[END_REF]. For the release case, the most widely used equation is the analytical solution of the Fick's second law for solute diffusion from a sphere to a well-stirred solution of finite volume [START_REF] Crank | The Mathematics of Diffusion[END_REF],

C t C eq = 1 - ∞ n=1 6α(1 + α) 9 + 9α + q 2 n α 2 exp - Dq 2 n R 2 t , (1)
where C t is the solute concentration in the bulk at time t, which is assumed to be uniform due to stirring. C eq is the expected concentration at equilibrium, i.e. as t → ∞. R is the radius of the sphere, D is the solute diffusion coefficient in the sphere, and α = V b /V is the ratio of the volumes of the bulk and the sphere, when assuming the partition coefficient is unity. q n is the n th non-zero positive root of tan q n = 3q n /(3 + αq 2 n ). C eq depends on the initial concentration C 0 within the sphere. It can be deduced from the mass conservation C 0 V = C eq (V b + V ) that leads to C eq = C 0 /(α + 1). For the absorption case, namely for diffusion from a well-stirred bulk solution to the capsule, the analytical solution is rather [START_REF] Crank | The Mathematics of Diffusion[END_REF],

C t C 0 = α 1 + α 1 + ∞ n=1 6 (1 + α) 9 + 9α + q 2 n α 2 exp - Dq 2 n t R 2 , (2)
where C 0 is, here, the initial concentration in the bulk. However, Eqs. ( 1) and ( 2) are not derived originally for a capsule with a composite core-shell structure, but simply for a homogeneous sphere (see Ref. [START_REF] Bielinski | Impact of the flow on mass transfer from particles: Biomedical applications[END_REF] for the mathematical derivation). The diffusion coefficient D is then considered by most authors as the effective diffusion coefficient of the overall capsule (core and membrane altogether). In this way, the fit enables the determination of an effective diffusion coefficient for the whole capsule, and not specifically for the membrane. However, the membrane permeability is a key parameter in designing, for example, capsules as bioreactors encapsulating cells and moderating the nutrients exchange rate via their membrane [START_REF] Ma | Core-shell hydrogel microcapsules for improved islets encapsulation[END_REF]. Only few authors have attempted to fit their experimental data with mathematical models derived for membranes. Kondo [START_REF] Kondo | Preparation and permeability characteristics of microcapsule membranes[END_REF] has used a model for a planar membrane by assuming negligible effects of the capsule curvature when the aspect ratio of the membrane thickness δ to the capsule radius R tends to zero (δ/R → 0). Henning et al. [START_REF] Henning | Characterizing permeability and stability of microcapsules for controlled drug delivery by dynamic NMR microscopy[END_REF] have measured concentration profiles using the NMR technique, and have fitted their experimental data with numerical solutions computed with the Matlab pdepe solver [24], instead of Eq. ( 2). The present article proposes a novel straightforward and correct methodology to characterise capsules permeability by an inverse approach, using solutions of the diffusion equation computed with the finite difference method. This article is organised as follows. The proposed method is explained in details in Sec. II, where the governing equations, the numerical scheme and its validation are presented. The efficiency of the suggested method is then demonstrated in Sec. III by applying it to numerical data and available experimental data in the literature (e.g. Refs. [START_REF] Koyama | Evaluation of mass-transfer characteristics in alginate-membrane liquid-core capsules prepared using polyethylene glycol[END_REF][START_REF] Rolland | Physicochemical properties of aqueous core hydrogel capsules[END_REF]). Finally, conclusions are drawn in Sec. IV.

II. METHOD

Governing equations -Diffusion of chemical species in a medium composed of an inner core and an outer domain separated by a shell that moderates mass exchange between these two compartments is governed by the Fick's second law in absence of advection,

∂c ∂t = ∇ • (D∇c) , ( 3
)
where t is the time and c the local concentration. D is the position-dependent diffusion coefficient taking different values in the core (0 ≤ r ≤ R c ), in the shell (R c < r ≤ R) and in the surrounding bulk (R < r ≤ R b ),

D(r) =      D c , 0 ≤ r ≤ R c D m , R c < r ≤ R D b , R < r ≤ R b , (4)
with r being the distance from the centre of the capsule. As the shell is a cross-linked gel or a porous rigid wall, its diffusivity to the solute D m is expected to be lower than the diffusivities in the core and the bulk, if these are liquids. We refer to the diffusion coefficient D m as the permeability of the capsule membrane to solute, see for example [START_REF] Kondo | Preparation and permeability characteristics of microcapsule membranes[END_REF][START_REF] Henning | Characterizing permeability and stability of microcapsules for controlled drug delivery by dynamic NMR microscopy[END_REF], while having in mind other definitions such as the ratio of the diffusion coefficient D m to the shell thickness δ [START_REF] Bielinski | Unsteady mass transfer from a core-shell cylinder in crossflow[END_REF]. All these definitions are adapted to cases where the shell thickness is not negligible compared to the size of the capsule δ/R ≫ 0.1. The problem presents a spherical symmetry. Thus, the diffusion is expected to be purely radial with c = c(r, t). Equation ( 3) is consequently expressed only as a function of its radial derivative terms,

∂c ∂t = 1 r 2 ∂ ∂r r 2 D ∂c ∂r = 2D r ∂c ∂r + ∂D ∂r ∂c ∂r + D ∂ 2 c ∂r 2 . ( 5
)
Initial condition -Diffusion from and to capsules is considered in the present study. For the release case, i.e. diffusion from the capsule to the bulk, the core of the capsule is initially loaded with a uniform concentration C 0 , and the concentration is zero elsewhere,

c(r, t = 0) = C 0 , 0 ≤ r ≤ R c 0, r > R c . (6)
This establishes a concentration gradient that triggers diffusion from the core to the bulk. For the absorption case, i.e. diffusion into the capsule, the core and the membrane are initially free of solute while the concentration is uniform and equal to C 0 in the bulk,

c(r, t = 0) = 0, r ≤ R C 0 , r > R . (7)
This condition, on the contrary, induces diffusion from the bulk to the capsule. Boundary conditions -The resulting spatial distribution of the concentration depends also on the boundary conditions. Here, unsteady continuous concentration and mass flux boundary conditions emerge in the absence of interfacial mass resistance at the core-shell interface

r = R c , c(R - c , t) = c(R + c , t) and D c ∂c ∂r r=R - c = D m ∂c ∂r r=R + c , (8)
and at the shell

-bulk interface r = R, c(R -, t) = c(R + , t) and D m ∂c ∂r r=R - = D b ∂c ∂r r=R + , (9)
where the superscripts (-) and (+) refer to the inner and outer sides of the interfaces, respectively. These boundary conditions relax the constraints of imposing perfect sink or well-stirred boundary conditions. These latter could, in principle, be fulfilled in an experiment, but not necessarily under operational conditions. They limit the computation to the capsule by excluding the bulk. The perfect sink conditions can be modeled with a Dirichlet boundary condition at the capsule surface, c(R, t) = 0, [START_REF] Fischer | Triggered release of loads from microcapsule-in-microcapsule hydrogel microcarriers: en-route to an "artificial pancreas[END_REF] and the well-stirred solution with the boundary condition,

∂c ∂t = -D m 3R 2 R 3 b -R 3 ∂c ∂r r=R , (11)
to be set at the surface of the capsule (see Ref. [START_REF] Bielinski | Impact of the flow on mass transfer from particles: Biomedical applications[END_REF] for the mathematical derivation of this equation). When the study considers three layers (core, shell and bulk), no mass flux is imposed at the domain edge (r = R b ),

∂c ∂r r=R b = 0. ( 12
)
The condition,

∂c ∂r r=0 = 0, ( 13
)
is set at the origin (r = 0) in all the simulations to fulfill the radial symmetry.

Numerical FD solver -The governing equations are solved numerically using the finite difference method (FD), while considering continuity of both the concentration and the mass flux at the inner and the outer interfaces of the shell. The discretisation of the equations is adopted from the scheme proposed by Hickson et al. [START_REF] Hickson | Finite difference schemes for multilayer diffusion[END_REF] for diffusion in one-dimensional multilayer slabs. It is adapted, here, to one-dimensional spherical coordinates that introduce additional terms. The diffusion equation in spherical coordinates, Eq. ( 5), is discretised at the inner points of each subdomain (represented as filled circles in Fig. 2) by a central finite difference approximation in space.

D i D i+1 j -3 j -2 j +2 j +3 j -1 j +1 j i i +1 ∆r !"#$% !"#$% r ∂Ω i FIG. 2.
Scheme showing the spatial discretisation in the vicinity of the interface ∂Ωi between two adjacent layers i and i + 1, where the solute diffusivities are respectively Di and Di+1. The inner points are represented with filled circles, and the point at the interface is shown as an empty circle. The space step ∆r is uniform in the whole domain.

The capsule is considered to be constituted of homogeneous layers, i.e. layers of constant diffusivity, which leads to,

∂c j ∂t = D i ∆r 2 j -1 j c j-1 -2c j + j + 1 j c j+1 , (14)
for the discretisation at the inner points of layer i. Here, c j is the local concentration at a given inner point j of the layer i, D i is the diffusion coefficient in the layer i, and ∆r the space step, which is set uniform in the whole domain. The discretisation of the domain proposed by Hickson et al. assumes one grid-point j lies on the interface ∂Ω i between the layers i and i + 1, as illustrated in Fig. 2. Taking a central difference scheme in space to discretise Eq. ( 5) at the interface ∂Ω i , and then taking first-order forward and backward differences for the spatial derivatives of the concentration gives,

∂c j ∂t = 1 2∆r 2 D i j -1 j 2 c j-1 -D i j -1 j 2 + D i+1 j + 1 j 2 c j + D i+1 j + 1 j 2 c j+1 . ( 15
)
Equations ( 14) and ( 15) can be written together under the form of a single linear system,

d dt                . . . c j-3 c j-2 c j-1 c j c j+1 c j+2 c j+3 . . .                =                    . . . . . . • • • j-3 j-2 χ i -2χ i j-1 j-2 χ i 0 0 0 0 • • • • • • 0 j-2 j-1 χ i -2χ i j j-1 χ i 0 0 0 • • • • • • 0 0 a -(a + b) b 0 0 • • • • • • 0 0 0 j j+1 χ i+1 -2χ i+1 j+2 j+1 χ i+1 0 • • • • • • 0 0 0 0 j+1 j+2 χ i+1 -2χ i+1 j+3 j+2 χ i+1 • • • . . . . . .                                   . . . c j-3 c j-2 c j-1 c j c j+1 c j+2 c j+3 . . .                , ( 16
)
where c j is the concentration at the interface ∂Ω i located at the node j (see Fig. 2),

χ i = Di ∆r 2 , χ i+1 = Di+1 ∆r 2 , a = χi 2 j-1 j 2
, and b = χi+1 2 j+1 j

2

. The concentration at iteration n + 1 is computed using a standard Euler scheme,

C n+1 = C n + ∆tMC n , ( 17
)
where C n is the column vector containing the values of the concentration at each grid-point at iteration n, ∆t is the timestep, and M the tridiagonal matrix in Eq. [START_REF] Leick | Alginate/poly-L-lysine capsules: mechanical properties and drug release characteristics[END_REF]. Here, M includes extra terms compared to the matrix given in Ref. [START_REF] Hickson | Finite difference schemes for multilayer diffusion[END_REF] because of the spherical coordinates for which is derived.

Validation -The FD scheme is validated by comparing the computed concentration profiles with semi-analytical solutions of Ref. [START_REF] Kaoui | Mechanistic modelling of drug release from multi-layer capsules[END_REF], with R c = 1.5 mm, δ = 0.2 mm, R b = 30 mm, D m = 0.5 × 10 -10 m 2 /s, and D c = D b = 3 × 10 -10 m 2 /s. The concentration profiles at times t = 5, 50 and 150 min are presented in Fig. 3a. For clarity reasons, the figure is truncated at r/R = 3 beyond which the concentration is almost zero. The membrane location is depicted by the light grey rectangle. The expected physics is well recovered with continuous and discontinuous variations of respectively the concentration and its derivative across the membrane. Perfect agreement with the semi-analytical solution is achieved, as well. The profiles are computed using the Matlab script diffusion FD.m provided in Ref. [START_REF] Bielinski | diffusion FD[END_REF]. For this simulation, the space step is set to ∆r = 10 -5 m which corresponds to 3000 elements. A grid independence study has been conducted to make sure the obtained results are reliable and invariant of grid resolution. Figure 3b 5) obtained with the proposed finite difference scheme (FD), while the green dashed lines represent the semianalytical solution taken from Ref. [START_REF] Kaoui | Mechanistic modelling of drug release from multi-layer capsules[END_REF]. The dark grey solid line gives the local solute diffusion coefficient D(r). The light grey rectangle depicts the location of the membrane. Excellent agreement is achieved between the FD and the semi-analytical solutions. (b) Concentration computed at the center of the capsule at t = 50 min for various grid resolutions. The result converges to a unique value at high grid resolutions.

shows the concentration at the capsule center (r = 0) at time t = 50 min computed for different grid resolutions using the same simulation parameters as in Fig. 3a. The results are not affected when refining the discretization further than ∆r = 10 -5 m (3000 elements), which is a good compromise between accuracy and computing time.

Permeability characterisation -The proposed method is based on an inverse analysis approach. The required inputs are the experimental data providing either the release or the uptake kinetics of a solute, and the geometrical characteristics of the studied capsule. A wide interval of the fitting parameter space (D m , D c ) is scanned, and the resulting output release or uptake curves are then compared against the experimental data to quantify the root mean square error (RMSE),

RMSE = 1 N N i=1 (C exp (i) -C num (i)) 2 , ( 18
)
where N is the number of the input experimental data points C exp , and C num the numerical solution. Only cases with D m ≤ D c are considered because diffusion in the membrane is expected to be slower than in the core, since the role of the membrane is to moderate mass exchange. The optimal values of the diffusivities (D m , D c ) are those leading to the minimal error between the numerical solution and the experimental input data. They are computed using the Matlab script compute permeability.m provided in Ref. [START_REF] Bielinski | compute permeability[END_REF].

III. RESULTS

A. Application based on numerical data

First, the method is tested on numerical data computed with a fully three-dimensional lattice Boltzmann simulation (LB), for which the diffusion coefficients are set and known in advance by the authors. Further details on the LB method can be found in Refs. [START_REF] Sukop | Lattice Boltzmann modeling[END_REF][START_REF] Krüger | The lattice Boltzmann method -Principles and practice[END_REF]. The geometrical parameters used in the LB simulation are R c = 1 mm, δ = 0.5 mm, and R b = 12 mm. The diffusion coefficients are D c = D b = 2 × 10 -10 m 2 /s and D m = 0.6 × 10 -10 m 2 /s corresponding to the diffusion of insulin in an aqueous solution and a polymeric membrane, respectively [START_REF] Morvan | Unsteady diffusion mass transfer in a microencapsulated islet of Langerhans for a bioartificial pancreas[END_REF]. Continuous boundary conditions, Eqs. ( 8) and ( 9), are considered. The resulting concentration field in the capsule and in its surrounding medium is given in Fig. 4a at various times. The release kinetics is quantified by 1 -Mt M0 , where M t = 4π R r=0 r 2 c(r, t)dr is the mass of the solute that remains within the capsule at time t. It is denoted by grey-colored square symbols in Fig. 4b. The characterisation method consists in taking these data (as if they were experimental data C exp ) and fit them with FD numerical solutions using the solute diffusivity in both the membrane and the core (D m , D c ) as fitting parameters, while holding constant the values of the geometrical parameters R c , δ and R b . Figure 4c reports the RMSE measuring how much the obtained FD solutions deviate from the LB solution. The red solid line gives the value of the membrane diffusivity D m which minimises the RMSE at a given core diffusivity D c . The global minimum of the RMSE is obtained at D c = 2.1 × 10 -10 m 2 /s and D m = 0.6 × 10 -10 m 2 /s (symbol ⊗). This estimated value of the membrane diffusivity is exactly equal to the value set in the LB simulation, and the small difference obtained for the core diffusivity represents a relative error of only 5%, which is satisfactory. The line D m = D c corresponds to the case of having a homogeneous sphere, as usually assumed in the literature. The local minimum of RMSE on this line gives the effective diffusivity D for the whole capsule. It is located at

D m = D c = D = 1.1 × 10 -10 m 2 /s (symbol ⊕).
The RMSE is strongly sensitive to changes in the membrane diffusivity D m , and it is weakly affected by variations of the core diffusivity D c , especially at low membrane diffusivities. This means that the release process is mainly controlled by the membrane permeability.

The release curve computed for D c = 2.1 × 10 -10 m 2 /s and D m = 0.6 × 10 -10 m 2 /s fits perfectly the input LB release data (see Fig. 4b), in contrast to the one of a homogeneous sphere. This example has allowed us to validate the proposed method and to demonstrate its efficiency to characterise the membrane permeability of capsules.

B. Application to experimental data of glucose release

Here, the proposed method is used to determine the membrane permeability to glucose of hydrogel alginate capsules, whose release kinetics has been measured experimentally by Rolland et al. [START_REF] Rolland | Physicochemical properties of aqueous core hydrogel capsules[END_REF]. In this study, the size of the capsules is R = 1.73 mm, and the thickness of the membrane varies within the range 45 µm < δ < 100 µm that leads to 0.026 < δ/R < 0.06. Only capsules of average membrane thickness δ = 50 µm have been used. They have been suspended in an aqueous solution whose glucose concentration is kept uniform thanks to stirring. For this case, the computational domain is limited only to two subdomain layers: the core and the shell excluding the bulk, where the concentration C t is forced to be uniform but yet varies over time. The well-stirred solution condition is modeled by Eq. [START_REF] Dembczynski | Characterization of small molecules diffusion in hydrogel-membrane liquid-core capsules[END_REF]. The initial concentration is set to C 0 = 0.1 g/mL in the whole capsule, i.e. in both the core and the membrane.

The experimental data of Rolland et al. are reported in Fig. 5a as square symbols. They show the evolution in time of the glucose concentration in the liquid surrounding the capsule. The fit of these experimental data, while scanning the values of the diffusion coefficients in both the core and the membrane, leads to the RMSE reported in value is about two times larger than the one in an infinitely dilute solution, that is 6.75 × 10 -10 m 2 /s. They attributed their observed increase in D to the stirring and motion of the capsule that speeds up the release. Here, the obtained value of D m is close to 6.75 × 10 -10 m 2 /s, which is expectedly small. The flow around the capsule induced by stirring is expected to enhance the release rate due to convection [START_REF] Bielinski | Unsteady mass transfer from a core-shell cylinder in crossflow[END_REF][START_REF] Kaoui | Flow and mass transfer around a core-shell reservoir[END_REF][START_REF] Wang | Enhancement of mass transfer from particles by local shear-rate and correlations with application to drug dissolution[END_REF], but not to significantly increase the intrinsic permeability of the capsule membrane.

C. Application to experimental data of glucose absorption

The present characterisation method can also be used in the case of solute absorption. Here, it is applied to glucose absorption by capsules made of alginate and polyethylene glycol using the experimental data of Koyama and Seki [START_REF] Koyama | Evaluation of mass-transfer characteristics in alginate-membrane liquid-core capsules prepared using polyethylene glycol[END_REF]. In this study, the inner core radius of the capsules is R c = 2.17 mm and the membrane thickness is δ = 160 µm. The prepared capsules have been immersed into a well-stirred solution of glucose that diffuses to the inner core of the capsules through their membrane. The authors have measured the variation of the glucose concentration in the surrounding bulk over time, which is reported in Fig. 6a as square symbols. The RMSE between these experimental data points and the numerical solution obtained by the FD method while varying the fitting parameters D c and D m is shown in Fig. 6b. As in the case of solute release, the RMSE is very sensitive to the variations of the membrane diffusivity D m . However, it is slightly affected by changes in the core diffusivity D c , especially at low membrane diffusivities. The absorption process is thus mainly controlled by the membrane permeability

IV. CONCLUSIONS

The proposed method has been applied to determine the permeability of capsules in the cases of release and absorption of a solute, while considering continuous and well-stirred boundary conditions at the capsule membrane. It determines the diffusion coefficient in both the membrane and the core, in contrast to the classically used approach that can only give an effective diffusion coefficient for the whole capsule. The computed diffusion coefficients in the membrane are smaller than the effective diffusion coefficients determined by the classical method. This is expected since the membrane is a porous medium whose role is to control and slow down the mass transfer. Furthermore, the present method recovers also the effective diffusion coefficients estimated by the classical method. Its ease-of-use makes the proposed method an efficient tool for determining accurately the capsule membrane permeability for drug delivery applications. It can also be used as a computer-aided tool in designing future capsules with desired release or absorption kinetics. The corresponding Matlab script is provided in Ref. [START_REF] Bielinski | compute permeability[END_REF].

  FIG.1. Scheme showing the structure of a core-shell capsule with radius R. Rc is the radius of the core and R b the radius of the bulk where the capsule is immersed. The shell is supposed to be homogeneous with a uniform thickness δ. The scheme also shows the mass transport directions of the release and absorption of a solute.

FIG. 3 .

 3 FIG. 3. (a)Concentration profiles computed at times t = 5, 50 and 150 min. The red solid line curves represent the numerical solution of Eq. (5) obtained with the proposed finite difference scheme (FD), while the green dashed lines represent the semianalytical solution taken from Ref.[START_REF] Kaoui | Mechanistic modelling of drug release from multi-layer capsules[END_REF]. The dark grey solid line gives the local solute diffusion coefficient D(r). The light grey rectangle depicts the location of the membrane. Excellent agreement is achieved between the FD and the semi-analytical solutions. (b) Concentration computed at the center of the capsule at t = 50 min for various grid resolutions. The result converges to a unique value at high grid resolutions.

FIG. 4 .

 4 FIG. 4. (a) Snapshots of the concentration field in three dimensions computed at various times by the lattice Boltzmann method (LB) in the case of insulin release from a capsule. The black solid lines delimit the core and the shell of the capsule. The domain is truncated at r/R b = 1/2 for clarity reasons. The indicated times are rounded down to the nearest hour. (b)The release kinetics curve. Squares correspond to the release kinetics computed by the LB method and used as input data for the characterisation. The green solid curve corresponds to the best fit of the LB data when considering the capsule with a composite structure, and the red solid curve corresponds to the best fit when assuming the capsule as a homogeneous sphere. (c) The RMSE between the LB and the FD solutions over a wide range of the fitting parameters Dc and Dm. The red solid line gives the value of the membrane diffusivity Dm that minimises the RMSE at a given core diffusivity Dc. The absolute minimum is obtained at Dc = 2.1×10 -10 m 2 /s and Dm = 0.6×10 -10 m 2 /s (symbol ⊗), and that match the values set in the LB simulation. These fit better the data than the solution obtained when considering the capsule as a homogeneous sphere with D = Dm = Dc = 1.1 × 10 -10 m 2 /s (symbol ⊕). The core radius, the membrane thickness and the bulk radius are Rc = 1 mm, δ = 0.5 mm and R b = 12 mm, respectively.

Fig. 5b .

 5b The minimum of the RMSE located on the line D m = D c (symbol ⊕) corresponds to the assumption that the capsule is a homogeneous sphere, as done by Rolland et al., whereas the minimum located in the area D m < D c (symbol ⊗) considers the composite core-shell nature of the capsule. The estimated value for a homogeneous sphere is D = 13.32 × 10 -10 m 2 /s, which is close to D = 14 × 10 -10 m 2 /s measured by Rolland et al. When taking into account the real composite structure of the capsule, the estimated diffusivities are rather D c = 17.11 × 10 -10 m 2 /s and D m = 7.54 × 10 -10 m 2 /s. The corresponding release curve fits the experimental data better than the one obtained for the homogeneous sphere, as demonstrated in Fig. 5a. Moreover, Rolland et al. have pointed out that their measured

FIG. 5 .FIG. 6 .

 56 FIG. 5. (a) The release kinetics of glucose from a hydrogel alginate capsule of size R = 1.73 mm and membrane thickness δ = 50 µm. Square symbols correspond to experimental data taken from Ref. [18], whereas the solid lines correspond to fits obtained with the present proposed method. (b) The RMSE error between the FD solution and the experimental data while varying the solute diffusion coefficient in both the core and the membrane. The global minimum of RMSE (⊗) gives the diffusion coefficients Dc = 17.11 × 10 -10 m 2 /s and Dm = 7.54 × 10 -10 m 2 /s that fit better the experimental data than the effective diffusion coefficient D = 13.32 × 10 -10 m 2 /s (⊕) obtained when assuming the capsule as a homogeneous sphere.

  . The global minimum (symbol ⊗) is obtained at D c = 76.40 × 10 -10 m 2 /s and D m = 2.96 × 10 -10 m 2 /s. The local minimum along the line D m = D c (symbol ⊕) that corresponds to a homogeneous sphere is located at D m = D c = D = 7.98 × 10 -10 m 2 /s, and it is very close to D = 7.9 × 10 -10 m 2 /s estimated by Koyama and Seki. The computed absorption curves are plotted in Fig. 6a. The curve of the capsule model with an internal composite structure again fits perfectly the experimental data.

ACKNOWLEDGEMENTS

The authors thank the Ministère de l'Enseignement Supérieur, de la Recherche et de l'Innovation (MESRI) and the Biomechanics and Bioengineering Laboratory (BMBI) for financial support.

DISCLOSURE STATEMENT

The authors report no conflict of interest.