Abstract
This paper faces two of the main drawbacks in networked control systems: bandwidth constraints and timevarying delays. The bandwidth limitations are solved by using multirate control techniques. The resultant multirate controller must ensure closed-loop stability in the presence of time-varying delays. Some stability conditions and a state feedback controller design are formulated in terms of linear matrix inequalities. The theoretical proposal is validated in two different experimental environments: a crane-based test-bed over Ethernet, and a maglev based platform over Profibus.
1. Introduction
Networked control systems (NCSs) [1, 2] are becoming a powerful research area because of introducing considerable advantages [3] (wiring reduction, easier and cheaper maintenance, cost optimization) in several kinds of control applications (teleoperation, supervisory control, avionics, chemical plants, etc). Nevertheless, as a consequence of sharing the same communication medium among different devices (sensor, actuator, controller), some problems such as time-varying delays [4–6], bandwidth limitations [7, 8], packet dropouts [9–13], packet disorder [14], and lack of synchronization [15] can arise in this kind of systems. So, the analysis and design of NCS is a complex problem and, usually, some simplifying assumptions are made.
In this work, neither packet dropouts nor packet disorder will be considered (as later detailed). In addition, devices will be assumed to be synchronized (by means of a suitable initial synchronization procedure or by implementing time-stamping techniques). Thus, only bandwidth limitations and time-varying delays will be faced. To solve these issues, some previous authors’ developments such as those in [4, 16–19] are revised and conveniently adapted and gathered together in the present work.
Regarding bandwidth constraints, this could be the case if the network configuration imposes a limitation in the control frequency (say, because of an excessive number of devices sharing the communication link). In this context, the consideration of a dual-rate controller [20–22] might be useful in terms of achievable performance [23, 24], since the controller can work at a faster rate than the network one which provides the measurements (a multirate input control (MRIC) structure will be considered, which is able to generate control actions for each sampled output).
As time-varying delays are assumed in an NCS, the control problem becomes a linear time-varying (LTV) one. Then, stability and design for arbitrary time-varying network delays must be carried out. Depending on what kind of information about the delay is provided, different stability analysis can be performed. So, if the probability density function of the network delay is unknown, a robust stability analysis must be proved. However, if the probability function is provided, stochastic stability can be analyzed. In this work, in order to prove any of these situations in a dual-rate NCS, linear matrix inequalities (LMIs) [25] will be considered. So, both LMIs for the robust case and for the probabilistic one (with extension to a multiobjective analysis) will be formulated. With respect to design approaches, a state feedback controller enunciated in terms of LMIs will be presented. The reader is referred, for example, to [6, 13] to find other LMI-based state feedback controller approaches. As well in terms of LMIs, in [11, 12] controllers are enunciated, and in [21] a multirate controller is proposed.
As a summary, the main novelties introduced by this work can be lumped as the followoing. (i)NCS analysis improvements: in [4, 16] control system stability is studied via LMIs. In the former work, a robust analysis is treated, whereas in the latter work, a stochastic analysis is presented. Both works use as a controller an output feedback one. In the present work, notation used in both approaches is unified and extended to the state feedback controller case. In addition, a hierarchical structure for the controller is contemplated. (ii)NCS design improvements: firstly, in [19], a multirate controller design procedure is introduced by splitting the controller into two sides—the slow-rate side and the fast-rate side. No shared communication medium is considered between both sides. Now, in the present work, from the previous idea a distributed controller for the proposed multirate NCS can be implemented (details in Section 2). Secondly, in [18], a single-rate state feedback controller to deal with time-varying delays is proposed. Now, this controller is adapted to be included in the slow-rate side of the multirate NCS. In addition, in order to reach a null steady-state error, an integral action can be added to this controller. The obtained results show a better control system performance than that achieved in [17].
The paper is divided into five parts: in Section 2 the dual-rate NCS scenario is exposed. Closed-loop realizations from lifted plant and controller are formulated in Section 3. By means of LMIs, for this kind of LTV systems, Section 4 proposes several stability analysis, and Section 5 introduces a state feedback controller. Finally, Section 6 presents stability and design results obtained for the test-bed platforms, and Section 7 enumerates the main conclusions.
2. Problem Scenario
Depending on the network configuration, three main options arise when integrating a dual-rate controller in an NCS. (i)The dual-rate controller is located at a remote side (with no direct link to the plant), and its fast-rate control actions can be sent from this side to the local actuator (directly connected to the plant) following a packet-based approach [26]. (ii)The dual-rate controller is located at the local side, being directly connected to the actuator [16]. (iii)The dual-rate controller is split into two subcontrollers [19] (a slow-rate one located at the remote side, and a fast-rate one situated at the local side); then a rude slow-rate control action is sent from the remote side to the local one in order to be refined and converted to a fast-rate control signal by the local subcontroller [4].
From the last option, Figure 1 represents a typical chronogram for this kind of dual-rate NCS. Both time-driven and event-driven policies are required [27]. The meaning of the encircled numbers is now detailed.(1)The sensor works in a time-triggered operation mode, sampling the process output at period (the measurement period). The output is sent through the network.(2)After a certain processing and propagation time has elapsed , the remote controller receives the packet.(3)Then, at the remote controller an event is triggered. As a consequence, after a computation delay , a slow-rate control action is generated and sent to the local controller, which is directly connected to the actuator.(4)After a certain processing and propagation time has elapsed , the local controller receives the packet.(5)Then, at the local side an event is triggered. Its main consequence is the application of N faster-rate control actions to the process. Such actions are scheduled to be applied taking into account the total delay: where, in this work, , being in order to avoid causality and packet order issues, and no sample losses. Note that the lumped delay in (2.1) can be adopted if the controller is a static one. In addition, if the total delay can be measured and finally compensated at the local side (as assumed in this work), (2.1) can also be considered in the control strategy. In any other case, it is inappropriate or impossible to lump all delays as one. Please, see for example in [28, 29] on how to deal with delays on different communication links separately.
Regarding the control actions, the first of them will be applied at the time of arrival of the packet ( time units after the measurement was taken). The remaining control actions, as they are not influenced by the network delay, will be applied every T time units (the control period), triggered by a fast-rate clock signal. Note, if the delay fulfills the first control actions will never be applied to the process.
3. Preliminaries and Notation
Consider a continuous linear time-invariant plant , which admits a state-space realisation , with suitable dimensions, fulfilling the followig: Being an arbitrary number, is denoted as It is well known [30, 31] that, in the case the input changes every time units, being constant in the intersample period (zero-order hold, ZOH) and the output is sampled synchronously to that input, the sampled output verifies the discrete-time equations:
When the input change and output sampling do not follow such a conventional sampling pattern, but they follow an arbitrary but periodic one with period , the discretization is a periodic linear time-varying discrete system. However, the process can be equivalently represented by a multivariable linear time-invariant discrete system with period and the so-called “lifted” input and output vectors, which are formed by stacking all the input and output signals; this methodology is denoted as “lifting” [32].
3.1. Lifted Plant Realization
Consider the above system (3.1) (as most physical system have , it will be assumed on the sequel) being subject to inputs at time , under ZOH conditions (i.e., input at time is held constant until time reaches , and , (see Figure 1). It is easy to show [31, 33] that the state at an arbitrary time is given by: and, from the above formula, a lifted realization can be computed when inputs are applied at times and outputs are read at times inside a metaperiod . Indeed, the discrete state equations comes from replacing in (3.4), and the output equation comes from replacing in (3.4) and multiplying by . In the following, as the above equations will be evaluated every seconds, notation will describe the input at time , and similarly, will denote the sample at time .
In a networked control framework, since is the last controller output from the previous sampling period (the controller will be assumed to apply (with delay ) a set of control actions (see Figure 1), an additional set of states must be added [18, 34]. So, incorporating the “memory” equation , and replacing by , (3.4) yields the following: The output equations would be built by stacking for all needed the following expression:
As described, an MRIC strategy is considered in this work, and hence only the first sampled output is needed to be sent to the controller. Then, for the sake of simplicity, let us describe the lifted plant model in this way: where , , , and being . Standard complete and matrices can be reviewed, for instance, in [32].
3.2. Controller and Closed-Loop Realization
In this paper, two different structures for the controller will be taken into account: a one-degree-of-freedom linear regulator , and a hierarchical controller .
For the first case, the regulator , two alternative cases can be treated as follows:(i) an output feedback regulator, whose lifted discrete realization will be [32]: where set-points are considered to be zero, so (being the loop error), (ii) a state feedback regulator, with a gain :
From the previous plant representation , its closed-loop connection to the output feedback regulator implies a dynamical system governed by this expression [31]:
For the state feedback regulator, the closed-loop realization will yield the following [31]:
Regarding the second regulator, the hierarchical one , it will be decomposed into two parts at different rates (remember Section 2). Different alternatives could be used to design each part. For brevity, let us consider and formulate the option developed in the second example of Section 6, where (i)fast-rate local subcontrollers are designed by means of robust control techniques [42], (ii)a coordinating, slow-rate remote subcontroller is designed using a state feedback approach (to be detailed in Section 5).
Then, the representation of each local subcontroller will be similar to (3.10). And, the remote subcontroller will yield an expression like in (3.11) but with these variations: where is the slow-rate control action, and is the resultant controller gain when considering the augmented state , which includes the overall local side state (controller + plant, remember (3.12)).
Let us denote the lifted expression , where, respectively, , , , are obtained like , , , in (3.8), but now considering the overall local side state (details omitted for brevity). Then, the closed-loop realization for the hierarchical control structure will be similar to (3.13), but now
4. Stability Analysis
As commented, time-varying delays can appear in an NCS framework. Thus, a variation in the instants where the outputs are measured () or those in which the input commands are presented to the plant () is expected. Let us denote the set of parameters that might vary from metaperiod to metaperiod as . Since matrices in (3.7) depend on , then , or , can vary from metaperiod to metaperiod. Further time variance occurs if the controller is also intentionally dependent on all or some of the parameters included in ; this is the case, for example, in a gain-scheduling approach [16, 35]. Subsequently, the closed-loop realization , , does depend on the above parameters, and then it will be replaced by in (3.12), by in (3.13), and by in (3.15), representing a discrete LTV system. For the sake of simplicity, only one of the closed-loop realizations () will be considered on the sequel. In addition, in this work the actual time-varying parameter used will be the network delay , that is, the delay of the first control action .
Three different scenarios will be studied as follows. (i)Consideration of arbitrary delay changes with unknown probability: a robust stability analysis will be needed. (ii)A probability density function of the network delay for each network situation is assumed known: a stochastic stability analysis can be independently carried out for each situation. (iii)Several network states are considered, which are defined by different probability density functions and different performance objectives: a multiobjective analysis will be developed.
4.1. Robust Analysis
In order to prove robust stability of the discrete LTV system: with a geometric decay rate (the geometric decay rate is a performance measure for nonlinear and LTV systems which guarantees that there exists so . When particularized to a discrete linear time-invariant system, the decay rate is the modulus of the dominant pole), , a common Lyapunov function must be found [18, 25, 36] so that (obviously, implies stability, given by the decrescence condition ). Replacing the closed-loop equations in (4.2), the Lyapunov decrescence condition can be written as the following LMI: where is a dummy parameter ranging in a set , where the time-varying parameters are assumed to take values in, and matrix is composed of decision variables to be found by the semidefinite programming solver. In this work, will be an interval (as defined in (2.1)).
If is an affine function of and is polytopic, then (4.3) can be checked with a finite number of LMIs. Otherwise, for bounded a dense enough gridding must be set up in order to approximately check for the above conditions. This procedure is denoted as LMI gridding [18, 36].
4.2. Probabilistic Analysis
Now, a probabilistic model of the network delay will be considered, so a probability density function is assumed known.
As the network delay is supposed to vary in a random way, stability of the closed-loop system will be analyzed in the mean square sense [37] (a particular case of the Martingales convergence theorem [38]), by means of the quadratic Lyapunov function (4.2), which will be shown to decrease in average. So, denoting as the statistical expectation, will tend to zero, and hence the state will converge to zero with probability one. The average descent to prove will be expressed as or, considering an average decay rate , the descent expression yields
Replacing the closed-loop equations in (4.2), the Lyapunov decrescence condition (4.5) can be written as the following probabilistic LMI: where , , and were defined after (4.3).
For a generic probability distribution, working with the above integral may be cumbersome. For bounded , a dense enough gridding in must be set up in order to approximately check for the above conditions. This procedure extends the LMI gridding in [18, 36] to a probabilistic case. Choosing a set of equally spaced values , so that , ( is an interval ), (4.6) can be approximately rewritten as which is a standard LMI to be solved by widely known methods [25, 39].
Note that the above results are more relaxed than those in the robust case. Indeed, in a probabilistic case there is only one LMI constraint (average decay) instead of one for each possible sampling period. In this way, temporal random increases of the Lyapunov function are tolerated as long as the average over time is decreasing. Hence, better results in stability analysis can be obtained; however, the gridding approach leaves intermediate points out of the analysis so, in rigor, the results are not valid unless the grid is very fine.
4.3. Multiobjective Analysis
The above idea can be extended to considering several possible network states, say , with different performance objectives. Each network state will be described by a probability density of the network delay , and a performance objective , . For example, two objectives can be considered: the first one can be defined by the probability function for an unloaded network (with a probability distribution around a “short” delay mean), and the other objective for a saturated network case (with a larger delay mean).
Then, the Lyapunov decrescence conditions can be written as the following probabilistic LMI (expressed, for computation, in its discrete approximation):
It is well known that the optimal result of multiobjective analysis will be a Pareto front with the optimal performance for a particular , being the rest fixed. If the performance bounds on the rest of constraints are made more restrictive, the resulting decreases. The reader is referred to [40] for basic ideas on multiobjective optimization.
Note that the use of the shared Lyapunov function in (4.8) proves stability of the networked control system for any probabilistic mixture of the considered network states (i.e., a network state whose probability density can be expressed as a convex combination of those in LMIs (4.8)).
5. State Feedback Controller Design
Apart from stability and decay rate, well-known LMI conditions can be set up for pole region placement, and norms, and so forth. The reader is referred to [25, 41] for details. In this section, a state feedback synthesis approach will be presented. The resultant controller can be used as a slow-rate subcontroller in a dual-rate framework.
As previously shown, from a lifted model , the control synthesis problem can be cast as a state-feedback one (3.14). But, as a consequence of time-varying network-induced delays, a different must be designed for each possible delay value (leading to, e.g., gain scheduling approaches [16, 35]). Another possibility is to achieve a unique stabilizing controller subject to any possible time-varying delay. In this case, an LMI gridding procedure can be considered. So, from [18], if there exist matrices and so that is verified for any , the feedback controller stabilizes with decay rate (which is the continuous-time equivalent exponential decay rate, that is, ), and is the associated Lyapunov function.
Due to the existence of time-varying delays, the lifted state could not be known when updating the control gain . This problem can be solved by adding any kind of state observer to the control strategy, for example, a nonstationary Kalman filter [18].
6. Experimental Results
6.1. A Crane-Based Platform over Ethernet: A Stability Analysis Example
In this section, a test-bed Ethernet environment is used to implement a dual-rate NCS, where the controller will be split into two parts. The proposed NCS includes the following devices (see Figure 2):(i) an industrial crane platform (to be controlled) equipped with three cc motors (to actuate each axis: ) and five encoders (to sense the three axis and two different angles). The motors are controlled by an analog signal in the range 1 V. The encoders provide a position measurement of 1 V/m. In this application only the -axis is actuated and sensed, whose behavior is modeled by
Details on the crane characteristics can be obtained at http://www.inteco.com.pl (3D crane apparatus). (ii) A local computer which is connected to the platform by means of a DAQ board, and where the local subcontroller is implemented. (iii) Two PLCs and one computer working as interference nodes in order to introduce different load scenarios. (iv) A switch shared by the previous devices to connect them to Ethernet. (v) A remote laptop computer where the remote subcontroller is implemented.
In this example, the controller will be a dual-rate PID one. Its parameters will be retuned according to Ethernet network delays, leading to a gain-scheduling proposal. In this case, the scheduling follows a Taylor-series-based approach (see [16] for details). An LMI analysis will be required to observe the stability benefits of the scheduled controller compared with the nominal unscheduled one. Finally, from experimental implementation, time response for both controllers will be obtained in order to observe transient behavior.
First of all, several experiments are carried out, where the number and complexity of the tasks developed by the interference nodes is modified in order to obtain different load scenarios ranging between the two extreme histograms on Figure 3. According to this figure, the output sampling time is chosen to be s, since the largest delay obtained is 0.39 s (delay bound: s), and then the requirement commented after (2.1), that is , is fulfilled.
Since the crane model (6.1) includes an integrator, a dual-rate PD is designed in order to achieve an overshoot of 8% and a settling time of 0.65 s (with no steady-state error). The resultant controller’s gains are: , , , (derivative filter).
6.1.1. Robust Stability Analysis
Let us suppose that no information about probability distribution is known. Then, the worst-case behavior of the proposed PD regulator can be assessed by means of the LMIs in (4.3). Testing different maximum delay bounds, the consequent results appear on Table 1.
As a conclusion, the proposed gain scheduled regulator improves worst-case performance for small delays (up to 0.2 s). In large delays, the approach used for retuning the PD parameters (based on Taylor-series) loses precision and results are similar (marginally worse) than those of a nonscheduled regulator. In fact, there are delay distributions involving delays larger than 0.3 s which might render the system unstable.
6.1.2. Probabilistic Stability Analysis: Extension to the Multiobjective Case
Now, information about probability distribution provided by experimental tests is taken into account. So, stability of the setup in probabilistic time-varying delays can be assessed. The LMI gridding in (4.7) can be carried out computing the closed-loop realization for the delay bound s . According to Figure 3, the number of grid points for the probability density approximation is taken as .
Two cases are analyzed as follows:(i)firstly, considering each network situation separately (a different Lyapunov function for each load scenario), the LMI in (4.7) is applied to each situation to obtain the minimum for which a feasible solution exists, (ii)secondly, a multiobjective analysis is performed by considering a unique Lyapunov function for both network load scenarios.
The second case is more conservative but allows stability guarantees for mixtures and random switching between both scenarios. The two proposed cases are somehow extreme situations from which would happen in a practical situation. If each of the network behaviors is very likely to remain active for a dwell time significantly longer than the loop’s settling time, then assumptions in case 1 will be closer to reality. If arbitrary, fast, network load changes were expected, then case 1 would be too optimistic and the analysis in case 2 would be recommended.
Regarding the first analysis, results are presented in Table 2, both for the scheduled PD and for the nominal one. In conclusion, the less the network is loaded, the better worst-case performance can be guaranteed. In addition, the scheduled approach shows better behavior than the nominal one.
Now, the second (multiobjective) study is carried out. Figure 4 shows a Pareto front that summarizes the analysis, which is developed by setting the decay-rate of one objective and optimizing the other’s one. As depicted, the decay-rates obtained in the previous study for the unloaded network case (here, the first objective) can not be now achieved, despite considering the highest decay-rate for the second objective (), that is, the loaded network scenario. However, if , the previous decay-rates for the loaded network case can be achieved. Finally, the figure reveals the scheduled approach outperforms again the nominal one.
In summary, from the analysis of both Tables 1 and 2, two main conclusions arise:(i)if the probability of large delays is low, the use of probabilistic information indicates (as intuitively expected) that the gain scheduling approach used in this example (based on Taylor-series) seems a sensible practical procedure, because of improving average performance. (ii)if no likelihood of (transient) instability is required, then the network must be reconfigured so the maximum delay does not exceed 0.3 s, or the initial controller specifications must be changed (reducing gains to improve robustness).
6.1.3. Control System Time Response
Since the previous figures indicate only stability and decay rate, to complete the study the control system time response is obtained. So, other performance differences (such as overshoot) can be evaluated.
Figures 5 and 6 present one of the different experiments tested for each network situation. As observed in the LMI analysis, the scheduled PD controller points out a better behavior than the nominal PD controller, being a more marked trend when working in a loaded network. So, Figures 5 and 6 show that the scheduled controller reduces the overshoot at least a 10% and up to a 40%, and the settling time up to a 60%, with a 30% on the average.
6.2. A Maglev-Based Platform over Profibus: A State Feedback Design Example
In this example, the position of a triangular platform assembled by joining three maglevs is hierarchically controlled by means of a dual-rate controller over Profibus. The proposed NCS (see Figure 7) includes:(i)a levitated platform with an equilateral triangle shape where each maglev is located at each corner. The maglevs provide position information from an infrared sensor array in V. The control signal is provided to a power amplifier, being in V. For more information about these levitators see http://www.xdtech.com, model ML-EA, (ii)a National Instruments CompactRio 9074 acting as local subcontrollers, (iii)a desktop PC acting as a remote subcontroller, (iv)a Profibus-DP network configured to work with a bus rate of 187.5 kBits/s, and with asynchronous operation mode. This enables sending a remote control action every 20 ms.
In this example, a standalone, fast-rate local subcontroller is designed for each maglev by using robust control techniques [42]. The coordinating, slow-rate remote subcontroller is a state feedback one, designed by using the LMI-gridding techniques presented in Section 5. So, the controller is assured to be robust in the presence of time-varying network-induced delays. Experimental time responses will validate this aspect.
After carrying out several experimental tests, network-induced time delays are measured (see histogram in Table 3). The main conclusion is that the most repeated round-trip time delay corresponds to a 10 ms period, with eventual delays at 5 ms and 15 ms. So, to design the state feedback controller, the grid of delay values to be considered will be ms, and hence (5.1) will be actually a collection of 3 LMIs. As the least delay value is 5 ms, then ms. And according to the bus rate, ms, and hence .
Now, the linearized state space for a generic maglev is presented as follows: where is intensity on levitator’s electromagnetic circuit, is the system output (i.e., a measure of airgap between levitated load and magnet taken with infrared sensors), , are resistance and inductance of the electromagnetic circuit for levitator , is mass of the levitated body, and , , , are constants of the magnetic levitator .
From this representation, and using the robust control toolbox in Matlab, a fast-rate, local subcontroller for each single maglev is designed ( ms). All the controllers are very similar, and hence one of them is now presented as follows:
If the coupled global platform model is obtained (details omitted for brevity; more information in [17]), the slow-rate, remote state feedback subcontroller will be designed. First, denoting the state, input, and output vectors as where , are angles of rotation of the levitated platform around - and -axes, respectively, and being . Then, the linearized state space for the coupled platform yields where Being, respectively, , moments of inertia around - and -axes of the levitated body.
From the previous model, the resulting state feedback controller in (3.14) obtained for a decay rate in (5.1) has the next gain matrix of dimensions :
As the design strategy contemplated in this example is the same than that used in [17], the controllers obtained in both cases show negligible differences.
6.2.1. Control System Time Response
Once the dual-rate controller is designed, the next experiment is carried out in the proposed network scenario (adding a nonstationary Kalman filter as a state observer).
The experiment starts with the platform in equilibrium point, as shown in Figure 8. In this figure, the top graphic shows the position error (center of mass), the middle one shows the control signal applied to the maglev, and the bottom one shows the supervision signal generated by the remote subcontroller and sent through the network to the local subcontroller. For clarity, only one of the three control and supervision signals is plotted. At time s, some load (a coin of 2 euros, 8.5 g) is applied. After a transient, the system reaches a new, stable equilibrium point, but with some position error (experiment 1). Possible disturbances introduced by coupling between the three maglevs are compensated by the coordinating remote subcontroller.
Next, as position error is presented, a new remote subcontroller that includes accumulated error in system state is developed. So the controller is designed considering where were introduced before (3.15), and the position error in will be zero in steady state [31].
Following this reasoning, the system state vector is expanded by adding the accumulated error for each one of the three maglevs. According to this new plant model, the feedback state controller is recalculated via LMI gridding, obtaining a new state feedback gain with dimensions (details omitted for brevity). As shown in Figure 8 (experiment 2), despite time-varying delays this new remote subcontroller (with the integral action) can keep the platform stable even with load variations (coin), improving control system performance with respect to that obtained in experiment 1 (which is related to [17]).
7. Conclusions
In this paper, in order to face arbitrary time-varying delays in a dual-rate NCS framework, different stability conditions and a state feedback design approach are presented in terms of LMIs. Multirate control techniques are proposed to avoid bandwidth limitations.
Regarding the stability conditions, three scenarios are treated: the robust case, the probabilistic case, and its extension to the multiobjective case. With respect to the state feedback controller, it is designed to assure robust stability for any possible time delay measured for the considered network.
Experimental results from two different dual-rate NCS implementations (a crane system over Ethernet, and a maglev-based platform over Profibus) validates the applicability of these LMI-based dual-rate control techniques.
Acknowledgments
The authors A. Cuenca, R. Pizá, and J. Salt are grateful to the Spanish Ministry of Education research Grants DPI2011-28507-C02-01 and DPI2009-14744-C03-03, and Generalitat Valenciana Grant GV/2010/018. A. Sala is grateful to the financial support of Spanish Ministry of Economy research Grant DPI2011-27845-C02-01, and Generalitat Valenciana Grant PROMETEO/2008/088.