I D Landau

D I Patrascu

Rejection of unknown multiple narrow band disturbances -A direct adaptive control approach

Keywords: direct adaptive control, internal model principle, Youla-Kucera parametrization, adaptive disturbance rejection, multiple narrow band disturbances

This paper presents a methodology for the rejection of multiple narrow band disturbances using a direct adaptive feedback control approach and without requiring the use of an additional transducer. The direct adaptive control scheme is based on the internal model principle and the use of the Youla-Kucera parametrization of the controller which allows a direct updating of the internal model in the controller. The performance of the methodology is evaluated on the Bodson model for the case of simultaneous application of several unknown time varying sinusoidal disturbances.

I. INTRODUCTION

One of the basic problems in control is the attenuation (rejection) of unknown disturbances without measuring them. The common framework is the assumption that the disturbance is the result of a white noise or a Dirac impulse passed through the "model of the disturbance". While in general one can assume a certain structure for such "model of disturbance", its parameters are unknown and may be time varying. This will require to use an adaptive approach. To be more specific, the disturbances considered are multiple narrow band or sinusoidal disturbances. These disturbances belong in fact to the class of "finite band disturbances". Furthermore for robustness reasons the disturbances should be located in the frequency domain within the regions where the plant has enough gain (see explanation in section III).

Solutions for this problem, provided that an "image" of the disturbance can be obtained by using an additional transducer, have been proposed by the signal processing community and a number of applications are reported ( [START_REF] Elliott | Active noise control. Noise / News International[END_REF], [START_REF] Elliott | Performance of feedforward and feedback systems for active control[END_REF], [START_REF] Beranek | Noise and Vibration Control Engineering: Principles and Applications[END_REF], [START_REF] Fuller | Active Control of Vibration[END_REF]). However, these solutions (inspired by Widrow's technique for adaptive noise cancellation ( [START_REF] Widrow | Adaptive Signal Processing[END_REF])) ignore the possibilities offered by feedback control systems and require an additional transducer. The basic idea is that a "well located" transducer can provide a measurement, highly correlated with the unknown disturbance. This information is applied to the control input of the plant through an adaptive filter (in general a Finite Impulse Response -FIR) whose parameters are adapted such that the effect of the disturbance upon the output is minimized. The disadvantages of this approach are:

• It requires the use of an additional transducer.

• Difficult choice for the location of this transducer (it is probably the main disadvantage). • It requires the adaptation of many parameters.

To achieve the rejection of the disturbance (at least asymptotically) without measuring it, a feedback solution can be considered. As mentioned earlier, the common framework is the assumption that the disturbance is the result of a white noise or a Dirac impulse passed through the "model of the disturbance" 1 . The present paper will focus on the case where the plant model is known but the disturbance model is unknown, since this is the situation encountered in many applications. Among the various approaches considered for solving this problem, the following ones may be mentioned:

1) Use of the internal model principle ( [START_REF] Francis | The internal model principle of control theory[END_REF], [START_REF] Johnson | Theory of disturbance-accomodating controllers[END_REF], [START_REF] Bengtsson | Output regulation and internal models -a frequency domain approach[END_REF], [START_REF] Tsypkin | Stochastic discrete systems with internal models[END_REF], [START_REF] Valentinotti | Adaptive Rejection of Unstable Disturbances: Application to a Fed-Batch Fermentation[END_REF], [START_REF] Ben Amara | Adaptive sinusoidal disturbance rejection in linear discrete-time systems -Part I: Theory[END_REF], [START_REF] Ben Amara | Adaptive sinusoidal disturbance rejection in linear discrete-time systems -Part II: Experiments[END_REF], [START_REF] Gouraud | Design of robust and frequency adaptive controllers for harmonic disturbance rejection in a single-phase power network[END_REF], [START_REF] Hillerstrom | Rejection of periodic disturbances with unknown period -a frequency domain approach[END_REF], [START_REF] Landau | Adaptive narrow band disturbance rejection applied to an active suspension -an internal model principle approach[END_REF]). 2) Use of an observer for the disturbance ( [START_REF] Marino | Robust adaptive compensation of biased sinusoidal disturbances with unknown frequency[END_REF], [START_REF] Ding | Global stabilization and disturbance suppression of a class of nonlinear systems with uncertain internal model[END_REF]).

3) Use of the "phase-locked" loop structure considered in communication systems ( [START_REF] Bodson | Adaptive algorithms for the rejection of sinusosidal disturbances with unknown frequency[END_REF], [START_REF] Bodson | Rejection of periodic distrubances of unknown and timevarying frequency[END_REF]). Of course, since the parameters of the disturbance model are unknown, all these approaches lead to an adaptive implementation which can be of direct or indirect type.

From the user point of view and taking into account the type of operation of existing adaptive disturbance compensation systems one has to consider two modes of operation of the adaptive schemes:

• Self-tuning operation (the adaptation procedure starts either on demand or when the performance is unsatisfactory and the current controller is frozen during the estimation/computation of the new controller parameters). • Adaptive operation (the adaptation is performed continuously and the controller is updated at each sampling). This paper explores the use of the internal model principle for the rejection of unknown time-varying multiple narrow band disturbances with the aim to obtain a direct adaptive control scheme.

Using the internal model principle, the controller should incorporate the model of the disturbance ( [START_REF] Francis | The internal model principle of control theory[END_REF], [START_REF] Johnson | Theory of disturbance-accomodating controllers[END_REF], [START_REF] Bengtsson | Output regulation and internal models -a frequency domain approach[END_REF], [START_REF] Tsypkin | Stochastic discrete systems with internal models[END_REF]). Therefore the rejection of unknown disturbances raises the problem of adapting the internal model of the controller and its re-design in real-time.

One way for solving this problem is to try to estimate in real time the model of the disturbance and re-compute the controller, which will incorporate the estimated model of the disturbance (as a pre-specified element of the controller). While the disturbance is unknown and its model needs to be

q B / A -d 1 / S 0 u(t) + + - y(t) Q A R 0 q B -d - -+ N / D p p p (t) 1 w(t) d(t)

Plant

Model Model

Adaptation Algorithm Ĉontroller Fig. 1. Direct adaptive control scheme for rejection of unknown disturbances estimated, one assumes that the model of the plant is known (obtained for example by identification). The estimation of the disturbance model can be done by using standard parameter estimation algorithms (see for example [START_REF] Landau | Techniques de modélisation récursive pour l'analyse spectrale paramétrique adaptative[END_REF], [START_REF] Ljung | System Identification -Theory for the User[END_REF]). This will lead to an indirect adaptive control scheme. The time consuming part of this approach is the redesign of the controller at each sampling time. The reason is that in many applications the plant model can be of very high dimension and despite that this model is constant, one has to re-compute the controller because a new internal model should be considered. This approach has been investigated in [START_REF] Bodson | Adaptive algorithms for the rejection of sinusosidal disturbances with unknown frequency[END_REF], [START_REF] Gouraud | Design of robust and frequency adaptive controllers for harmonic disturbance rejection in a single-phase power network[END_REF], [START_REF] Hillerstrom | Rejection of periodic disturbances with unknown period -a frequency domain approach[END_REF]. However, by considering the Youla-Kucera parametrization of the controller (known also as the Q-parametrization), it is possible to insert and adjust the internal model in the controller by adjusting the parameters of the Q polynomial (see figure 1). It comes out that in the presence of unknown disturbances it is possible to build a direct adaptive control scheme where the parameters of the Q polynomial are directly adapted in order to have the desired internal model without recomputing the controller (polynomials R 0 and S 0 in figure 1 remain unchanged). The number of the controller parameters to be directly adapted is roughly equal to the number of parameters of the denominator of the disturbance model. In other words, the size of the adaptation algorithm will depend upon the complexity of the disturbance model.

This paper focuses on the direct feedback adaptive control for the case of unknown and time-varying multiple narrow band disturbances. The direct adaptive control scheme to be presented( [START_REF] Landau | Adaptive narrow band disturbance rejection applied to an active suspension -an internal model principle approach[END_REF]) takes advantage of the Youla-Kucera parametrization for the computation of the controller. A related approach has been considered in [START_REF] Valentinotti | Adaptive Rejection of Unstable Disturbances: Application to a Fed-Batch Fermentation[END_REF] for an application to a chemical reactor but a theoretical analysis of the scheme is not provided. A related paper is also [START_REF] Ben Amara | Adaptive sinusoidal disturbance rejection in linear discrete-time systems -Part I: Theory[END_REF], [START_REF] Ben Amara | Adaptive sinusoidal disturbance rejection in linear discrete-time systems -Part II: Experiments[END_REF] where the application field is the active noise control in an acoustic duct. However the case of rejection of unknown multiple narrow band disturbances has not been discussed in any of the above mentioned references.

The paper is organized as follows. Section II is dedicated to a brief review of the plant, disturbance and controller representation as well as of the Internal Model Principle. Some robustness issues are addressed in section III. The direct adaptive control schemes for disturbance rejection is presented in section IV. Simulation results are presented in section V. Some concluding remarks are given in section VI.

II. PLANT REPRESENTATION AND CONTROLLER

STRUCTURE

The structure of a linear time invariant discrete time model of the plant (used for controller design) is:

G(z -1 ) = z -d B(z -1 ) A(z -1 ) = z -d-1 B * (z -1 ) A(z -1 ) , (1)
with: d = the plant pure time delay in number of sampling periods ;

A = 1 + a 1 z -1 + • • • + a n A z -n A ; B = b 1 z -1 + • • • + b n B z -n B = q -1 B * ; B * = b 1 + • • • + b n B z -n B +1 ,
where A(z -1 ), B(z -1 ), B * (z -1 ) are polynomials in the complex variable z -1 and n A , n B and n B -1 represent their orders 2 . The model of the plant may be obtained by system identification. Details on system identification of the plant models considered in this paper can be found in [START_REF] Landau | Digital Control Systems -Design, Identification and Implementation[END_REF], [START_REF] Landau | Direct controller order reduction by identification in closed loop[END_REF].

Since in this paper we are focused on regulation, the controller to be designed is a RS-type polynomial controller ( [START_REF] Landau | Adaptive control[END_REF], [START_REF] Landau | Digital Control Systems -Design, Identification and Implementation[END_REF]) -see also figure 1.

The output of the plant y(t) and the input u(t) may be written as:

y(t) = q -d B(q -1 ) A(q -1 ) • u(t) + p 1 (t) ; (2)
S(q -1 ) • u(t) = -R(q -1 ) • y(t) , (3)
where q -1 is the delay (shift) operator (x(t) = q -1 x(t + 1)) and p 1 (t) is the resulting additive disturbance on the output of the system. R(z -1 ) and S(z -1 ) are polynomials in z -1 having the orders n R and n S , respectively, with the following expressions:

R(z -1 ) = r 0 + r 1 z -1 + . . . + r n R z -n R = R ′ (z -1 ) • H R (z -1 ) (4) S(z -1 ) = 1 + s 1 z -1 + . . . + s n S z -n S = S ′ (z -1 ) • H S (z -1 ) (5)
where H R and H S are pre-specified parts of the controller (used for example to incorporate the internal model of a disturbance or to open the loop at certain frequencies).

We define the following sensitivity functions:

• Output sensitivity function (the transfer function between the disturbance p 1 (t) and the output of the system y(t)):

S yp (z -1 ) = A(z -1 )S(z -1 ) P(z -1 ) ; (6)
• Input sensitivity function (the transfer function between the disturbance p 1 (t) and the input of the system u(t)):

S up (z -1 ) = - A(z -1 )R(z -1 ) P(z -1 ) , (7)
where

P(z -1 ) = A(z -1 )S(z -1 ) + z -d B(z -1 )R(z -1 ) = A(z -1 )S ′ (z -1 ) • H S (z -1 ) +z -d B(z -1 )R ′ (z -1 ) • H R (z -1 ) (8)
defines the poles of the closed loop. In pole placement design, P(z -1 ) is the polynomial specifying the desired closed loop poles and the controller polynomials R(z -1 ) and S(z -1 ) are minimal degree solutions of [START_REF] Ding | Global stabilization and disturbance suppression of a class of nonlinear systems with uncertain internal model[END_REF] where the degrees of P, R and S are given by:

n P ≤ n A + n B + d -1, n S = n B + d -1 and n R = n A -1.
Using the equations ( 2) and

(3), one can write the output of the system as:

y(t) = A(q -1 )S(q -1 ) P(q -1 ) • p 1 (t) = S yp (q -1 ) • p 1 (t) . ( 9
)
For more details on RS-type controllers and sensitivity functions see [START_REF] Landau | Digital Control Systems -Design, Identification and Implementation[END_REF].

Suppose that p 1 (t) is a deterministic disturbance, so it can be written as

p 1 (t) = N p (q -1 ) D p (q -1 ) • δ (t) , (10)
where δ (t) is a Dirac impulse and N p (z -1 ), D p (z -1 ) are coprime polynomials in z -1 , of degrees n N p and n D p , respectively. In the case of stationary disturbances the roots of D p (z -1 ) are on the unit circle. The energy of the disturbance is essentially represented by D p . The contribution of the terms of N p is weak compared to the effect of D p , so one can neglect the effect of N p .

Internal Model Principle:

The effect of the disturbance given in [START_REF] Elliott | Performance of feedforward and feedback systems for active control[END_REF] upon the output:

y(t) =
A(q -1 )S(q -1 ) P(q -1 )

• N p (q -1 ) D p (q -1 ) • δ (t) , (11)
where D p (z -1 ) is a polynomial with roots on the unit circle and P(z -1 ) is an asymptotically stable polynomial, converges asymptotically towards zero if and only if the polynomial S(z -1 ) in the RS controller has the form:

S(z -1 ) = D p (z -1 )S ′ (z -1 ) . (12)
In other terms, the pre-specified part of S(z -1 ) should be chosen as H S (z -1 ) = D p (z -1 ) and the controller is computed using [START_REF] Ding | Global stabilization and disturbance suppression of a class of nonlinear systems with uncertain internal model[END_REF], where P, D p , A, B, H R and d are given 3 . Using the Youla-Kucera parametrization (Qparametrization) of all stable controllers ([3], [START_REF] Tsypkin | Stochastic discrete systems with internal models[END_REF]), the controller polynomials R(z -1 ) and S(z -1 ) get the form:

R(z -1 ) = R 0 (z -1 ) + A(z -1 )Q(z -1 ) ; (13) S(z -1 ) = S 0 (z -1 ) -z -d B(z -1 )Q(z -1 ) . (14)
The (central) controller (R 0 , S 0 ) can be computed by poles placement (but any other design technique can be used).

Given the plant model (A, B, d) and the desired closed-loop poles P one has to solve:

P(z -1 ) = A(z -1 )S 0 (z -1 ) + z -d B(z -1 )R 0 (z -1 ) . ( 15
)
3 Of course it is assumed that D p and B do not have common factors.

Equations ( 13) and ( 14) characterize the set of all stabilizable controllers assigning the closed loop poles as defined by P(z -1 ) (it can be verified by simple calculations that the poles of the closed loop remain unchanged). For the purpose of this paper Q(z -1 ) is considered to be a polynomial of the form:

Q(z -1 ) = q 0 + q 1 z -1 + . . . + q n Q z -n Q . ( 16
)
To compute Q(z -1 ) in order that the controller incorporates the internal model of the disturbance one has to solve the diophantine equation:

S ′ (z -1 )D p (z -1 ) + z -d B(z -1 )Q(z -1 ) = S 0 (z -1 ) , (17)
where D p (z -1 ), d, B(z -1 ) and S 0 (z -1 ) are known and S ′ (z -1 ) and Q(z -1 ) are unknown. Equation ( 17) has a unique solution for

S ′ (z -1 ) et Q(z -1 ) with: n S 0 ≤ n D p + n B + d -1, n S ′ = n B + d -1, n Q = n D p -1 .
One sees that the order n Q of the polynomial Q depends upon the structure of the disturbance model.

III. ROBUSTNESS CONSIDERATIONS

As it is well known, the introduction of the internal model for the perfect rejection of the disturbance (asymptotically) will have as effect to raise the maximum value of the modulus of the output sensitivity function S yp . This may lead to unacceptable values for the modulus and the delay margins if the controller design is not appropriately done (see [START_REF] Landau | Digital Control Systems -Design, Identification and Implementation[END_REF]). As a consequence, a robust control design should be considered assuming that the model of the disturbance is known, in order to be sure that for all situations an acceptable modulus margin and delay margin are obtained.

On the other hand at the frequencies where perfect rejection of the disturbance is achieved one has S yp (e -jω ) = 0 and S up (e -jω ) =

A(e -jω ) B(e -jω ) .

Equation (18) corresponds to the inverse of the gain of the system to be controlled. The implication of equation ( 18) is that cancellation (or in general an important attenuation) of disturbances on the output should be done only in frequency regions where the system gain is large enough. If the gain of the controlled system is too low, |S up | will be large at these frequencies. Therefore, the robustness vs additive plant model uncertainties will be reduced and the stress on the actuator will become important. Equation ( 18) also implies that serious problems will occur if B(z -1 ) has complex zeros close to the unit circle (stable or unstable zeros) at frequencies where an important attenuation of disturbances is required. It is mandatory to avoid attenuation of disturbances at these frequencies.

Since on one hand we would not like to react to very high frequency disturbances and on the other hand we would like to have a good robustness it is often wise to open the loop at 0.5 f s ( f s is the sampling frequency) by introducing a fixed part in the controller H R (q -1 ) = 1 + q -1 (for details see [START_REF] Landau | Digital Control Systems -Design, Identification and Implementation[END_REF] and section II).

IV. DIRECT ADAPTIVE CONTROL FOR DISTURBANCE

ATTENUATION

The objective is to find an estimation algorithm which will directly estimate the parameters of the internal model in the controller in the presence of an unknown disturbance (but of known structure) without modifying the closed loop poles. Clearly, the Q-parametrization is a potential option since modifications of the Q polynomial will not affect the closed loop poles. In order to build an estimation algorithm it is necessary to define an error equation which will reflect the difference between the optimal Q polynomial and its current value.

In [START_REF] Tsypkin | Stochastic discrete systems with internal models[END_REF], such an error equation is provided and it can be used for developing a direct adaptive control scheme. This idea has been used in [START_REF] Valentinotti | Adaptive Rejection of Unstable Disturbances: Application to a Fed-Batch Fermentation[END_REF], [START_REF] Ben Amara | Adaptive sinusoidal disturbance rejection in linear discrete-time systems -Part I: Theory[END_REF], [START_REF] Ben Amara | Adaptive sinusoidal disturbance rejection in linear discrete-time systems -Part II: Experiments[END_REF], [START_REF] Landau | Adaptive narrow band disturbance rejection applied to an active suspension -an internal model principle approach[END_REF]. Using the Q-parametrization, the output of the system in the presence of a disturbance can be expressed as: y(t) = A(q -1 )[S 0 (q -1 )-q -d B(q -1 )Q(q -1 )] P(q -1 )

• N p (q -1 ) D p (q -1 ) • δ (t) = S 0 (q -1 )-q -d B(q -1 )Q(q -1 ) P(q -1 ) • w(t) , (19)
where w(t) is given by (see also figure 1):

w(t) = A(q -1 )N p (q -1 ) D p (q -1 ) • δ (t) = A • y(t) -q -d • B • u(t) . ( 20
)
In the time domain, the internal model principle can be interpreted as finding Q such that asymptotically y(t) becomes zero. Assume that one has an estimation of Q(q -1 ) at instant t, denoted Q(t, q -1 ). Define ε 0 (t + 1) as the value of y(t + 1) obtained with Q(t, q -1 ). Using [START_REF] Landau | Techniques de modélisation récursive pour l'analyse spectrale paramétrique adaptative[END_REF] one gets:

ε 0 (t + 1) = S 0 (q -1 ) P(q -1 ) w(t + 1) -q -d B * (q -1 ) P(q -1 ) Q(t, q -1 )w(t) (21) One can define now the a posteriori error (using Q(t + 1, q -1 )) as:

ε(t + 1) = S 0 (q -1 ) P(q -1 ) w(t + 1) - q -d B * (q -1 ) P(q -1 ) Q(t + 1, q -1 ) • w(t) (22)
Replacing S 0 (q -1 ) from the last equation by [START_REF] Landau | Direct controller order reduction by identification in closed loop[END_REF] one obtains

ε(t + 1) = [Q(q -1 ) -Q(t + 1, q -1 )] q -d B * (q -1 ) P(q -1 ) w(t) + v(t + 1) (23)
where v(t) = S ′ (q -1 )D p (q -1 ) P(q -1 )

• w(t) = S ′ (q -1 )A(q -1 )N p (q -1 )

P(q -1 ) • δ (t)
is a signal which tends asymptotically towards zero. Define the estimated polynomial Q(t, q -1 ) as: Q(t, q -1 ) = q0 (t) + q1 (t)q -1 + . . . + qn Q (t)q -n Q and the associated estimated parameter vector : θ (t) = [ q0 (t) q1 (t) . . . qn Q (t)] T . Define the fixed parameter vector corresponding to the optimal value of the polynomial Q as: θ = [q 0 q 1 . . . q n Q ] T . Denote:

w 2 (t) = q -d B * (q -1 ) P(q -1 )
• w(t) [START_REF] Valentinotti | Adaptive Rejection of Unstable Disturbances: Application to a Fed-Batch Fermentation[END_REF] and define the following observation vector:

φ T (t) = [w 2 (t) w 2 (t -1) . . . w 2 (t -n Q )] . (25)
Equation ( 23) becomes

ε(t + 1) = [θ T -θ T (t + 1)] • φ (t) + v(t + 1) . ( 26
)
One can remark that ε(t) corresponds to an adaptation error ( [START_REF] Landau | Adaptive control[END_REF]).

From equation ( 21) one obtains the a priori adaptation error:

ε 0 (t + 1) = w 1 (t + 1) -θ T (t)φ (t) ,
with

w 1 (t + 1) = S 0 (q -1 ) P(q -1 ) • w(t + 1) ; ( 27
)
w 2 (t) = q -d B * (q -1 ) P(q -1 )
• w(t) ; ( 28)

w(t + 1) = A(q -1 ) • y(t + 1) -q -d B * (q -1 ) • u(t) , ( 29
)
where B(q -1 )u(t + 1) = B * (q -1 )u(t).

The a posteriori adaptation error is obtained from [START_REF] Marino | Robust adaptive compensation of biased sinusoidal disturbances with unknown frequency[END_REF]:

ε(t + 1) = w 1 (t + 1) -θ T (t + 1)φ (t) .
For the estimation of the parameters of Q(t, q -1 ) the following parameter adaptation algorithm is used ( [START_REF] Landau | Adaptive control[END_REF]):

θ (t + 1) = θ (t) + F(t)φ (t)ε(t + 1) ; (30)
ε(t + 1) = ε 0 (t + 1) 1 + φ T (t)F(t)φ (t) ; (31)
ε 0 (t + 1) = w 1 (t + 1) -θ T (t)φ (t) ; (32)
F(t + 1) = 1 λ 1 (t)   F(t) - F(t)φ (t)φ T (t)F(t) λ 1 (t) C) + φ T (t)F(t)φ (t)   . ( 33
) 1 ≥ λ 1 (t) > 0; 0 ≤ λ 2 (t) < 2 (34)
where λ 1 (t), λ 2 (t) allow to obtain various profiles for the evolution of the adaption gain F(t) (for details see [START_REF] Landau | Adaptive control[END_REF], [START_REF] Landau | Digital Control Systems -Design, Identification and Implementation[END_REF]).

In order to implement this methodology for disturbance rejection (see figure 1), it is supposed that the plant model

z -d B(z -1 ) A(z -1 )
is known (identified) and that it exists a controller [R 0 (z -1 ), S 0 (z -1 )] which verifies the desired specifications in the absence of the disturbance. One also supposes that the degree n Q of the polynomial

Q(z -1 ) is fixed, n Q = n D p -1,
i.e. the structure of the disturbance is known.

The following procedure is applied at each sampling time for adaptive operation:

1) Get the measured output y(t + 1) and the applied control u(t) to compute w(t + 1) using (29). 2) Compute w 1 (t + 1) and w 2 (t) using ( 27) and (28) with P given by [START_REF] Johnson | Theory of disturbance-accomodating controllers[END_REF]. 3) Estimate the Q-polynomial using the parametric adaptation algorithm (30) -(33). 4) Compute and apply the control (see figure 1): S 0 (q -1 )u(t + 1) = -R 0 (q -1 )y(t + 1) -Q(t, q -1 )w(t + 1) (35) For the self tuning operation of the adaptive scheme, the estimation of the Q-polynomial starts once the level of the output is over a defined threshold. A parameter adaptation algorithm (30)-(33) with decreasing adaption gain is used and the estimation is stopped when the adaption gain is below a pre-specified level 4 . During estimation of the new parameters, the controller is kept constant. The controller is updated once the estimation phase is finished. For a stability analysis of this scheme see [START_REF] Landau | Adaptive narrow band disturbance rejection applied to an active suspension -an internal model principle approach[END_REF].

V. A SIMULATED EXAMPLE

The rejection of unknown disturbance using the adaptive control scheme proposed in Section IV is illustrated for the case of a model taken from [START_REF] Bodson | Adaptive algorithms for the rejection of sinusosidal disturbances with unknown frequency[END_REF] [START_REF] Bodson | Adaptive algorithms for the rejection of sinusosidal disturbances with unknown frequency[END_REF] where ẏ(t) = -a 0 y(t)+ b 0 u(t)+ v(t).The transfer function is:

H(s) = 1 0.01s + 1 (36)
In our case the disturbance will be represented by a sum of two, respectively four sinusoids with time-varying frequency applied to the output of the plant , so we shall consider n Dp = 4 and n Dp = 8 that mean for Q -parametrization n Q = 3 respectively n Q = 7. For a realistic simulation we have considered a random noise signal -0.1 ≤ v(t) ≤ 0.1 applied on the output of system plant. The implementation of the controller is in discrete time. The selected sampling frequency is f s = 250Hz The nominal digital controller [R 0 , S 0 ] (without the internal model of the disturbance) has been designed using pole placement with calibration of sensitivity functions. The dominant poles for closed loop are a pair of complex poles fixed at the frequency 136.36 rad/sec 5 with a damping ξ = 0.8 The considered disturbances applied on the output of system have the form: The parameter adaptation algorithm used for rejection of those disturbances has a number of parameters according with the complexity of d 1 (t) and d 2 (t) (see above for details). For the case of two sinusoids we use the following simulation protocol:

d 1 (t) = 0.
• The system operates in close loop with the noise v(t) applied all the time(0 ≤ t ≤ 30). • For 1 ≤ t < 10 the disturbance d 1 (t) = 0.6 sin( f 11 •2π)+ 0.6 sin( f 12 • 2π) with f 11 = 10Hz and f 12 = 15Hz is applied (each sinusoid has magnitude 0.6). 4 The magnitude of the adaptation gain gives an indication upon the variance of the parameter estimation error -see for example [START_REF] Landau | Adaptive control[END_REF]. 5 The frequency 136. and f 12 = 15Hz is applied again (20 ≤ t ≤ 30). The estimation algorithm works permanently and the controller is recomputed at each sampling period. The method for parameter adaptation is forgetting variable factor switching to constant trace (for details see [START_REF] Landau | Adaptive control[END_REF], [START_REF] Landau | Digital Control Systems -Design, Identification and Implementation[END_REF]) and the degree of Q -polynomial is 3. The time domain responses in open loop and in closed loop are presented in figure 2. One can see that the adaptation transient is negligible when the disturbance characteristics are changing. The level of the residual output is close to the level of the noise. We can see that the algorithm works well for both sums of sinusoids.: The spectral densities of the output in open loop, respectively in closed loop using the adaptive rejection are presented in Fig. 3. We remark that the magnitudes of the two frequencies are attenuated by 30 dB bringing the residual disturbance almost at the level of the noise (in the absence of noise the attenuation is almost 40 dB). In the case of four sinusoids the simulation conditions are:

• The system operates in close loop with the noise v(t) applied all the time(0 ≤ t ≤ 30.

• For 1 ≤ t ≤ 30 the disturbance d 2 (t) with f 21 = 4Hz, f 22 = 9Hz, f 23 = 13Hz, f 24 = 17Hz is applied. The estimation algorithm works permanently and the controller is recomputed at each sampling period. The method for parameter adaptation is forgetting variable factor switching to constant trace and the degree of Q -polynomial is output in open loop, respectively in closed loop using the adaptative rejection are presented in Fig. 5. We remark that the amplitudes for all frequencies are attenuated by almost 30 dB bringing the residual disturbance close to the level of the noise (in the absence of noise the attenuation is more than 35 dB). The results obtained in a realistic simulation including the presence of a measurement noise are very good both in terms of attenuation and adaption transients.

The approach presented extent the methodology already used in practice for active vibration control in the presence of a single narrow band disturbance.

For future research it will be interesting:

• to evaluate comparatively this approach with the other approaches in terms of complexity and performance. • to include an anti windup procedure .

• to study the robustness of the scheme with respect to plant model uncertainties

  36 rad/sec (22 Hz) assures a step response in close loop almost equal to that in open loop • For 10 ≤ t < 20 a change in frequencies occurs. The new frequencies are: f 11 = 5Hz respectively f 12 = 8Hz. • At t = 20, the first sum of sinusoids with f 11 = 10Hz

Fig. 2 .

 2 Fig. 2. Time domain results of adaptive rejection for a sum of 2 sinusoids (upper figure: open loop, lower figure: closed loop).

Fig. 3 .

 3 Fig. 3. Power spectral density for a sum of 2 sinusoids.

7 .Fig. 4 .

 74 Fig. 4. Time domain results of adaptive rejection for a sum of 4 sinusoids(upper figure: open loop, lower figure: closed loop).

Fig. 5 .

 5 Fig. 5. Power spectral density for adaptive rejection of a sum of 4 sinusoids.

  11 , f 12 , f 21 , f 22 , f 23 , f 24 are the frequencies in Hz for different sinusoids.

	6 sin( f 11 2πt) + 0.6 sin( f 12 2πt)	(37)
	d 2 (t) = 0.6 sin( f 21 2πt) + 0.6 sin( f 22 2πt)
	+0.6 sin( f 23 2πt) + 0.6 sin( f 24 2πt)	(38)
	where the f

Throughout the paper it is assumed that the order of the disturbance model is known but the parameters of the model are unknown (the order can be estimated from data if necessary).

The complex variable z -1 will be used for characterizing the system's behavior in the frequency domain and the delay operator q -1 will be used for describing the system's behavior in the time domain.


Jannik Dreier
email: jannik.dreier@imag.fr

Florian Kerschbaum
email: florian.kerschbaum@sap.com

Practical Privacy-Preserving Multiparty Linear Programming Based on Problem Transformation

Cryptographic solutions to privacy-preserving multiparty linear programming are slow. This makes them unsuitable for many economically important applications, such as supply chain optimization, whose size exceeds their practically feasible input range. In this paper we present a privacy-preserving transformation that allows secure outsourcing of the linear program computation in an efficient manner. We evaluate security by quantifying the leakage about the input after the transformation and present implementation results. Using this transformation, we can mostly replace the costly cryptographic operations and securely solve problems several orders of magnitude larger.

I. INTRODUCTION

Linear Programming (LP) can be used to solve many practical optimization problems, e.g. supply chain master planning [START_REF] Stadler | Supply chain management and advanced planning[END_REF]. Many practical problems are distributed and require protection of the input data. For example, in supply chain master planning, the participating companies need to exchange information about production costs and capacities. This is very sensitive data, since it directly impacts the negotiation position. Consequently no master planning solution will be adopted in practice that reveals this data [START_REF] Pibernik | An approach to inter-domain master planning in supply chains[END_REF].

Secure multi-party computation (SMC) [START_REF] Ben-Or | Completeness theorems for non-cryptographic fault-tolerant distributed computation[END_REF][START_REF] Goldreich | How to play any mental game[END_REF][START_REF] Shamir | How to share a secret[END_REF][START_REF] Yao | Protocols for Secure Computations[END_REF]] offers a cryptographic solution to the problem. Several parties can jointly compute a function, e.g. LP, without disclosing anything except what can be inferred by a party's input and output. In theory this offers an ideal solution to the conflict posed by these problems. There even exist a number of specialised protocols for secure LP [START_REF] Catrina | Secure multiparty linear programming using fixed-point arithmetic[END_REF][START_REF] Li | Secure and private collaborative linear programming[END_REF][START_REF] Toft | Solving linear programs using multiparty computation[END_REF].

Nevertheless these solutions suffer from a prohibitively bad computational performance. We estimate that our prototypical implementation of Toft's protocol [START_REF] Toft | Solving linear programs using multiparty computation[END_REF] requires 7 years in order to solve our use case LP problem with 282 variables. Supply chain planning problems can easily reach 4 million variables [START_REF] Braun | Optimization with grid computing[END_REF] and these are only single company planning problems while we consider cross-organizational optimization. These size of problems are currently solved every day by non-secure LP solvers. Assuming an average computation complexity of O(n 3 ) for LP and that Moore's Law continues c 2011 IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media, including reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works. to hold for the time being, it would still take roughly 80 years until these problems could be solved securely as fast as they can be solved non-securely today. Clearly, if we want to solve these problems today, we need a different approach.

In [START_REF] Vaidya | Privacy-preserving linear programming[END_REF] Vaidya presents a different approach to the problem. Instead of implementing a LP solver using SMC the problem is randomly transformed and then the problem is solved using a non-secure LP solver. Unfortunately, as Bednarz already points out in [START_REF] Bednarz | Hiccups on the road to privacy-preserving linear programming[END_REF], Vaidya's transformation is incorrect and may lead to solutions which are not admitted in the original problem. In this paper we present a new, correct transformation. Furthermore we do not settle for an informal assessment of its security, but instead evaluate it in the framework of leakage quantification. Finally, we present performance results from our implementation.

Let L be an instance of a LP problem where parties X i have input x i , e.g. variables, constraints or costs. The approach of [START_REF] Catrina | Secure multiparty linear programming using fixed-point arithmetic[END_REF][START_REF] Li | Secure and private collaborative linear programming[END_REF][START_REF] Toft | Solving linear programs using multiparty computation[END_REF] is to implement a LP solver for L using SMC where x i is the input of each party. Let g be our privacy-preserving transformation, then L = g(L) is another instance of a LP problem, but L reveals little information about L. We can outsource the solution of L to an untrusted service provider, e.g. the cloud. Either a single party computes g on its data, e.g. a single company instance of supply chain planning, or we can securely compute g using SMC for multiple parties. We show that computing g is much more efficient than implementing a LP solver and even present an optimized secure computation protocol for it. The solution computed by the cloud on L then needs to be transformed by g to the solution of L. Figure 1 depicts our approach to privacy-preserving linear programming using problem transformation.

In summary this paper contributes

• a novel privacy-preserving transformation for linear programs regardless of their partitioning.

• a proof that differently from previous work the transformation is correct, i.e. an optimal solution to the transformed problem always corresponds to an optimal solution of the original problem. • an analysis that the transformation is secure in the framework of leakage quantification. In our use case example the chances of guessing e.g. the optimal values x are less than 2.30 • 10 -1409 . • a secure computation protocol to efficiently compute the transformation in a multiparty case. • a theoretical analysis and performance results of our implementation that show that the transformation is efficient. In our use case example computing the solution only required 25 minutes (compared to 7 years for Toft's protocol). The remainder of this paper is structured as follows. In the next Section we review related work. In Section III we present our main result, the privacy-preserving transformation. We prove the correctness of the transformation in Section IV and analyze its performance in Section V. In Section VI we summarize our results from leakage quantification of this transformation. We show the protocol for securely computing the transformation in Section VII. The results from our use case example are described in Section VIII. We conclude the paper in Section IX.

II. RELATED WORK

Our work is related to secure multi-party LP solvers [START_REF] Catrina | Secure multiparty linear programming using fixed-point arithmetic[END_REF][START_REF] Li | Secure and private collaborative linear programming[END_REF][START_REF] Toft | Solving linear programs using multiparty computation[END_REF], secure outsourcing [START_REF] Atallah | Securely outsourcing linear algebra computations[END_REF][START_REF] Atallah | Secure outsourcing of scientific computations[END_REF][START_REF] Bednarz | Hiccups on the road to privacy-preserving linear programming[END_REF][START_REF] Du | Privacy-preserving cooperative scientific computations[END_REF][START_REF] Mangasarian | Privacy-Preserving Horizontally-Partitioned Linear Programs[END_REF][START_REF] Mangasarian | Privacy-Preserving Linear Programming[END_REF][START_REF] Vaidya | Privacy-preserving linear programming[END_REF] and leakage quantification [START_REF] Braun | Quantitative notions of leakage for one-try attacks[END_REF][START_REF] Smith | On the Foundations of Quantitative Information Flow[END_REF][START_REF] Wibmer | Leakage Quantification of Cryptographic Protocols[END_REF].

A. Secure Multi-Party LP Solvers

As mentioned above, different distributed LP solvers based on secure multi-party computation have been developed. Toft [START_REF] Toft | Solving linear programs using multiparty computation[END_REF] used a distributed secure Simplex algorithm based on secret sharing. This algorithm is usable for any number of participants and provides information-theoretic security. Yet it turned out to be too slow in practice for realistically sized problems in many applications such as Supply Chain Management.

Li and Atallah [START_REF] Li | Secure and private collaborative linear programming[END_REF] also used a distributed Simplex algorithm. In contrary to Toft [START_REF] Toft | Solving linear programs using multiparty computation[END_REF] they concentrated on the case where the data is shared among only two parties. It is unclear if their algorithm can be efficiently extended to the multi-party case.

Catrina and de Hoogh [START_REF] Catrina | Secure multiparty linear programming using fixed-point arithmetic[END_REF] developed a more efficient version of the distributed Simplex algorithm using fixed-point arithmetics. However, their solution still is only suitable for small inputs as it suffers from a higher asymptotic and worst-case complexity than our transformation.

B. Secure Outsourcing

Vaidya [START_REF] Vaidya | Privacy-preserving linear programming[END_REF] used a transformation approach which is the basis for our transformation. He supposed a scenario where one party holds the objective function and the other party the constraints. As a consequence his transformation does not protect the entire Linear Program, which is not appropriate for most applications. Additionally, he did not provide a formal security analysis and his work suffers from correctness problems [START_REF] Bednarz | Hiccups on the road to privacy-preserving linear programming[END_REF].

Mangasarian [START_REF] Mangasarian | Privacy-Preserving Horizontally-Partitioned Linear Programs[END_REF][START_REF] Mangasarian | Privacy-Preserving Linear Programming[END_REF] proposed similar transformations for special input distributions (horizontal or vertical partitioning). Supply chain optimization problems have a more complex distribution which combines both cases. Additionally, his approach for horizontal partitioning is limited to equality constraints. This is not sufficient for supply chain optimization either, as we have to deal with capacity constraints expressed as inequality constraints. Our transformation addresses both shortcomings and is suitable for more general data distributions. Moreover he does not provide a formal security analysis.

Furthermore there is work on related problems such as matrix multiplication and linear systems of equations [START_REF] Atallah | Securely outsourcing linear algebra computations[END_REF][START_REF] Atallah | Secure outsourcing of scientific computations[END_REF][START_REF] Du | Privacy-preserving cooperative scientific computations[END_REF].

C. Leakage Quantification

To give a formal security assessment of our transformation, we use methods of Leakage Quantification which measure how much information about the secret input is revealed to an attacker. In particular, we concentrate on the metric of "Multiplicative Advantage" which was first developed -but not called so -by Smith [START_REF] Smith | On the Foundations of Quantitative Information Flow[END_REF] based on entropy loss using a special type of min-entropy. Later on, this was analyzed more closely by Braun et al. [START_REF] Braun | Quantitative notions of leakage for one-try attacks[END_REF] who proposed another, simpler definition -the one we use in this paper. Compared to other metrics such as entropy loss or channel capacity using the standard Shannon-entropy, it is particularly expressive for onetry attacks as it gives a worst-case evaluation. To assess complex operations we make use of some theorems by Wibmer et al. [START_REF] Wibmer | Leakage Quantification of Cryptographic Protocols[END_REF] on the leakage of combined channels.

III. THE TRANSFORMATION

Linear Programming is a standard tool in business optimization. A Linear Program (LP) consists of a set of unknown variables x, a linear target function c(x) representing the costs which shall be minimized (or equivalently the gain which has to be maximized) and a set of constraints (linear equalities or inequalities):

min c T x s.t. M 1 x = b 1 M 2 x ≤ b 2 x ≥ 0
We use a positive monomial matrix1 Q to hide c (as proposed by [START_REF] Vaidya | Privacy-preserving linear programming[END_REF][START_REF] Bednarz | Hiccups on the road to privacy-preserving linear programming[END_REF]):

min c T QQ -1 x M 1 QQ -1 x = b 1 M 2 QQ -1 x ≤ b 2 Q -1 x ≥ 0
and a positive vector r to hide x

min c T Q(Q -1 x + r) M 1 Q(Q -1 x + r) = b 1 + M 1 Qr M 2 Q(Q -1 x + r) ≤ b 2 + M 2 Qr (Q -1 x + r) ≥ r For z = Q -1
x + r and a strictly positive diagonal matrix S2 we have

min c T Qz M 1 Qz = b 1 + M 1 Qr M 2 Qz ≤ b 2 + M 2 Qr Sz ≥ Sr Then we define c T = c T Q, M =   M 1 Q 0 M 2 Q A -S   , b =   b 1 + M 1 Qr b 2 + M 2 Qr -Sr  
where A is a permutation matrix3 representing slackvariables 4 . This allows us to rewrite the program as follows:

min c T s z s s.t. M z s = b z s ≥ 0
where c s is c with added zeros for the slack-variables and z s is the variable vector (z with added slack-variables). To hide the contents of M and b we use a nonsingular matrix P and with M = P * M and b = P * b we have

min c T s z s.t. M z s = b z s ≥ 0
As z = Q -1 x + r, the resulting x can be obtained from z by calculating x = Q(z -r).

IV. CORRECTNESS

For the transformation to be useful in practice, we have to ensure that the transformed program can still be used to find a solution to the original problem. We show (for more details see the technical report [START_REF] Dreier | Practical Secure and Efficient Multiparty Linear Programming Based on Problem Transformation[END_REF]) that any optimal (i.e. with minimal cost) and feasible (with respect to the constraints) solution to the transformed problem corresponds to an optimal and feasible solution of the original problem using the following two lemmas. Lemma 1. A solution x is feasible in the original problem if and only if z = Q -1 x + r is a feasible solution to the transformed problem.

Proof: This is true by construction. The multiplication by P does not change the solution set as P is invertible. Q is monomial (which gives correctness for this part of the transformation as shown by Bednarz et al. [START_REF] Bednarz | Hiccups on the road to privacy-preserving linear programming[END_REF]) and r is positive to not interfere with z ≥ 0. As S is a strictly positive diagonal matrix, the multiplication by S is actually a multiplication of each row with a positive scalar which leaves the inequalities untouched.

Lemma 2. Let z be the solution that minimizes c T z in the transformed LP. Then x = Q(z -r) minimizes c T x in the original problem.

Proof: We will use a proof by contradiction. Suppose that there is a solution x with c T x < c T x. Then

c T x < c T Q(z -r) ⇔ c T x < c T Qz -c T Qr ⇔ c T QQ -1 x + c T Qr < c T Qz ⇔ c T (Q -1 x + r) < c T z (1)
Apparently z = Q -1 x + r is a valid solution for transformed LP with strictly lower cost than z. Hence z is not optimal.

V. PERFORMANCE

Performance is a key aspect of our approach: The transformation has to be significantly less complex than solving the problem itself, otherwise outsourcing would not represent an efficiency gain to the clients.

From a theoretical viewpoint this is easy to see, as the most expensive operation during the transformation is the multiplication by P . Let m 1 and m 2 denote the number of rows of M 1 and M 2 respectively, and n the number of variables. Then the multiplication using a naive algorithm is in O((

m 1 + m 2 + n) 2 × (2n + m 2 )
). This is more efficient than the simplex algorithm (on which the secure LP solvers [START_REF] Catrina | Secure multiparty linear programming using fixed-point arithmetic[END_REF][START_REF] Li | Secure and private collaborative linear programming[END_REF][START_REF] Toft | Solving linear programs using multiparty computation[END_REF] are built on) with exponential worst case complexity [START_REF] Klee | How good is the simplex algorithm? Inequalities[END_REF]. Even in the average case or compared to interior point methods performance will be better due to smaller constants.

We discuss practical performance in our use case example in Section VIII.

VI. SECURITY

We analyze the security of our transformation in the following setting: An attacker wants to obtain the input data, i.e. the original LP. He knows the transformed LP and some abstract facts about the input, for example that the input values are within a certain range.

Ideally we would like to give a classical security proof. However, as our transformation is based on disguising using random noise, our transformation will probably leak some information. This makes the classical cryptological security definitions unsuitable. Nevertheless leaking some information can be acceptable for many applications as long as the leakage is small and bounded. In return, we get dramatic performance improvements which can be more important than perfect security. Nevertheless, we want to formally evaluate this leakage, and to show that it is "small and bounded". This is possible using Leakage Quantification methods based on information theory.

Due to space limitations, we are not able to discuss the entire analysis here (for details see our Technical Report [START_REF] Dreier | Practical Secure and Efficient Multiparty Linear Programming Based on Problem Transformation[END_REF]). We limit ourselves to defining the metric "Multiplicative Advantage" and presenting the results of our analysis, i.e. the bounds on the leakage of different parts of the transformation. In the last part of this section we analyze possible attacks.

A. Multiplicative Advantage

The disguising transformation is modeled as a communication channel as in information theory, and we try to measure by how much the knowledge of the output increases the chances of guessing the input.

Definition 1 (Channel). A discrete, noisy and memoryless channel C is given by

• A finite set X = x 1 , ..., x n called the input alphabet, • a finite set Y = y 1 , ..., y n called the output alphabet, • for each x ∈ X a random variable C|x that takes values in Y.
In this model the input x ∈ X represents the data to hide, the output y ∈ Y is the "encrypted" or "disguised" data that is passed into the cloud. The definition of the input alphabet X is important as it reflects the attackers "abstract knowledge" about the input, for example that the input is a value within a certain range.

Braun et al. [START_REF] Braun | Quantitative notions of leakage for one-try attacks[END_REF] defined the notion of "advantage" based on a similar notion ("vulnerability") by Smith [START_REF] Smith | On the Foundations of Quantitative Information Flow[END_REF]. For a random variable X characterizing the input distribution let

P R priori (X) = max i p(X = x i )
denote the a priori probability of a right guess (i.e. the probability of an attacker correctly guessing the input if it has not yet seen the output of the channel). Similarly let

P R posteriori (X) = j max i (p(y j |x i )p(X = x i ))
denote the a posteriori probability of a right guess (i.e. the probability of an attacker correctly guessing the input after having seen the output of the channel). Then we define

L (X) = P R posteriori (X) P R priori (X)
the factor by which the knowledge of the output of the channel increases the chances of a right guess for a given X5 .

Definition 2 (Multiplicative Advantage). The Multiplicative Advantage or Multiplicative Leakage of a channel C is defined as

L(C) = max X L (X)
L(C) is independent of the input distribution, which is convenient if this distribution is not known. Additionally L (X) has the advantage of attaining its maximum for a uniform input distribution [START_REF] Braun | Quantitative notions of leakage for one-try attacks[END_REF][START_REF] Wibmer | Leakage Quantification of Cryptographic Protocols[END_REF] which yields

L(C) = j max i p(y j |x i )
This result is very useful in practice as p(y j |x i ) are the entries of the channel matrix. Furthermore Wibmer et al. [START_REF] Wibmer | Leakage Quantification of Cryptographic Protocols[END_REF] showed that multiplicative leakage can be used to easily bound the leakage of compositions of independent channels which simplifies the analysis of complex protocols. In our technical report [START_REF] Dreier | Practical Secure and Efficient Multiparty Linear Programming Based on Problem Transformation[END_REF] we calculate bounds on the leakage of several basic operations (permutations, scalar and matrix multiplication) which we now use to analyze the overall leakage of our transformation.

B. Establishing Bounds on the Leakage

Using the bounds on the leakage of the basic operations, we establish bounds on the leakage of the three main parts of the transformation: The cost function, the variables and the constraints.

1) Leakage of the cost function c: In our transformation c is hidden by a monomial matrix Q, i.e. a permutation followed by a multiplication of each entry with a scalar value. Wibmer et al. [START_REF] Wibmer | Leakage Quantification of Cryptographic Protocols[END_REF] showed that the leakage of two independent channels arranged one after the other is bounded by the minimum of both leakages. Hence we obtain

L(c l n,k ) ≤ min n+l l+1 + 1 k , (A(n,l)+k)! k! * (A(n,l))!
where c l n,k denotes the channel for a k-dimensional input vector c hidden by a multiplication by a monomial matrix and A(n, l) is the number of different values of f = l i=1 f i with f i ∈ {1, . . . , n}. The first value is the leakage of a scalar multiplication, the second value bounds the leakage of the permutation (cf. [START_REF] Dreier | Practical Secure and Efficient Multiparty Linear Programming Based on Problem Transformation[END_REF]). Note that this is an evaluation of c only, although Q appears in other places in the transformation as well (in the constraint matrix, on the right hand side and in z). When evaluating all parts together, this results in partly dependent channels. However, this simplification is realistic as the other occurrences are well-hidden, as we discuss later on.

2) Leakage of the variables x: The variable vector z s consists of two parts: The non-slack variables z and the slack-variables. The non-slack variables contain values z = Q -1 x + r.

Wibmer et al. [START_REF] Wibmer | Leakage Quantification of Cryptographic Protocols[END_REF] showed that for a channel y = rx + r with r < r ≤ r m and x ≤ x m the leakage is bounded between

log(x m + 2) ≤ L("y = rx + r ") ≤ 2 log(x m + 2)
for r m → +∞. The channel "z = Q -1 x + r" can be modeled as a permutation followed by the "y = rx + r "-channel as Q resp. Q -1 is a monomial matrix and thus we have

L("z = Q -1 x + r") ≤ min 2 log(x m + 2), (x m + k)! k! * x m ! P M 1 Q, M 2 Q S r b 1 , b 2 w.h.p. M 1 , M 2 , c (permuted)

Fig. 2. Overview of a possible attack

As mentioned before, the matrix Q is also used in the hiding of c. This is not a problem, since the vector r prevents the search for common factors (see cryptanalysis below).

The slack variables contain potentially correlated values as well (S * Q -1 x), but since they are permuted, it is hidden which variable belongs to which constraint. This further complicates exploiting correlation by common factors.

3) Leakage of the Constraints (M 1 , M 2 , b 1 and b 2 ): The leakage of the constraints is difficult to estimate due to the number and complexity of related transformation steps. We can however argue that their security is essentially based on P , as the following analysis shows.

On the one hand, breaking the hiding with P almost completely reveals all other hiding operations (see Figure 2). If the attacker knows P , he can obtain S, which allows him to calculate r. This gives access to b 1 and b 2 . As Q is a monomial matrix, each row of M {1,2} Q is multiplied with the same factor which he can guess with very high probability when searching for common factors (see the cryptanalysis for details). This then leaks a permutation of M 1 , M 2 and c, which may be a complete break if the attacker knows something about the structure of M 1 or M 2 which allows him to undo the permutation. Hence we can say that the constraints are at most as safe as P .

On the other hand they are unlikely less safe than P , since conversely the knowledge of the constraints noticeably increases the chances of guessing P .

We now give an estimation of the leakage of P . To simplify the analysis, we use the following channel model: The input alphabet X are invertible square matrices of dimension (m 1 + m 2 + k) × (m 1 + m 2 + k) (where m 1 and m 2 are the number of rows of M 1 and M 2 ; k is the number of variables). The channel output is

P * B = P * B' 0 A
where A is a permutation matrix of dimension (m 2 + k) × (m 2 + k), and B a random matrix of dimension (m 1 + m 2 + k) × (k + 1). B represents the parts of M which contain M 1 Q, M 2 Q and -S as well as b (which can be seen as just another column of the matrix). The general idea is that A is responsible for most of the leaked information about P , and that attackers have very limited knowledge about the other parts of the matrix, so that modeling them as random is reasonable. This model allows us to use another theorem by Wibmer et al. [START_REF] Wibmer | Leakage Quantification of Cryptographic Protocols[END_REF] on parallel channels. We split the channel C P into a first channel C P,1 which calculates P * B , and a second channel C P,2 which calculates

P * 0 A
Then the theorem gives that max {L(C P,1 ), L(C P,2 )} ≤ L(C P ) ≤ L(C P,1 )L(C P,2 )

C P,1 corresponds to a matrix multiplication channel M k,n,l [START_REF] Dreier | Practical Secure and Efficient Multiparty Linear Programming Based on Problem Transformation[END_REF]. The leakage of C P,2 can be bounded as follows

L(C P,2 ) ≤ m2+k i=1 ((n + 1) m1+m2+l -i) (m 2 + k)!
when P contains values between 0 and n. This is formally shown in the technical report [START_REF] Dreier | Practical Secure and Efficient Multiparty Linear Programming Based on Problem Transformation[END_REF].

C. Cryptanalysis

In this Section we discuss several attacks on P and Q, the most important parts of the transformation.

1) Attacks on P : The most promising approaches to attack P are based on the structure of the matrices, in particular the slack-variables as the following example illustrates. Consider these constraints:

M * x ≤ b x ≥ 0
When adding slack-variables we obtain (I is the identity matrix)

M I * x = b x ≥ 0
If we multiply by P the result is

P * M P * x = P * b x ≥ 0
This means that the slack-variables completely reveal P .

Unfortunately the most general type of matrix we can use instead of the identity matrix I is a monomial matrix A as A and A -1 have to be positive [START_REF] Kaczorek | Positive 1D and 2D systems[END_REF]. Yet this is only partly helpful as the knowledge of P * A allows an attacker to calculate

M = (P * A) -1 * (P * M ) = A -1 * P -1 * P * M = A -1 * M As A is a monomial matrix, A -1 * M is a permutation of M
where each column is multiplied with the same value -which can be obtained with very high probability when searching for common factors. Not transforming the inequalities into equalities is not an option, as this would reduce our choice of P to monomial matrices for correctness reasons. This would turn P * M into an easy target for factorization attacks and thus not provide enough security.

To solve the problem, we treat equality constraints separately. This reduces the size of A, as we need less slackvariables. If A is smaller than P , it only leaks a permutation of some columns of P and not the entire matrix. However, this requires that in the input problem some of the constraints are equations, which is not always the case. Yet in our target application Supply Chain Management most of the constraints are actually equations, which turns this into a very well-suited solution.

If the input problem contains only inequality constraints, we can still improve security by assigning non-zero costs to slackvariables and permuting the variables to hide them among the real variables. A possible way to do this without changing the optimal solution is by analyzing the dual problem [START_REF] Sierksma | Linear and Integer Programming: Theory and Practice[END_REF].

2) Attacks on Q: Q is a positive monomial matrix, i.e. it can be written as Q = D * E where D is a diagonal matrix and E a permutation matrix. Thus

Q -1 = (D * E) -1 = E -1 * D -1 = E T * D -1 where D =    D 1 • • • 0 . . . . . . . . . 0 • • • D n    , D -1 =    1 D1 • • • 0 . . . . . . . . . 0 • • • 1 Dn   
Q appears in several places in the transformation. This can be an issue if the resulting products of Q and Q -1 with M 1 , M 2 , c etc. allow the search for common factors. Since two random numbers have an asymptotic probability 6/π 2 ≈ 61% of being coprime [START_REF] Kiltz | Secure computation of the mean and related statistics[END_REF], this search is very likely to be successful if we have access to several numbers containing the same factor. Thus we have to ensure that there are no two or more values which contain the same factor. We now analyze all occurrences of Q inside the transformation.

• The cost function c = c T * Q = c T * (D * E) = c 1 * D 1 c 2 * D 2 • • • c n * D n * E
contains all the factors of D, but each factor appears only once inside one of the entries. • M {1,2} Q has the following structure

M 1 * Q = M 1 * (D * E) =    M 1,1 * D 1 • • • M 1,n * D n . . . . . . M m1,1 * D 1 • • • M m1,n * D n    * E
This is open for attack as there are lots of values containing the same factor. However, since inside M there is -S below (cf. Eq. III), M = P * M has a different structure. For example

M 1,1 = m1 i=1 P 1,i * (M 1 Q) i,1 + m2 i=1 P 1,i+m1 * (M 2 Q) i,1 -P 1,m1+m2+1 * S 1
contains a common factor D j in all summands except for the last. This is enough to make the values unusable for factorization attacks. • M {1,2} Qr is secure against factorization attacks since

(M 1 Qr) 1 = n i=1 (M 1 Q) 1,i * r i where each (M 1 Q) 1,i contains a different D j .
• z = Q -1 x + r is protected against factorization attacks by r. • The slack-variables could also leak the factors inside Q as the i-th constraint concerning S has the form -S i * z i + Z si = -S i * r i which gives

Z si = S i * (z i -r i ) = S i * ((Q -1 x) i + r i -r i ) = S i * (Q -1 x) i
This would allow the search for common factors. However, the slack-variables are permuted independently of Q so that it is unknown which variable belongs to which constraint. To find possible matches an attacker could search for pairs of variable whose difference is smaller than the maximum of r. However, this is made harder by

S as Z si -z i = S i * z i -S i * r i -z i = (S i -1)z i -S i * r i (
and not -r i as it would be the case if S i = 1 for all i). Overall only in c the factors inside Q are easily accessible. Thus the search for common factors is difficult, if possible at all, which is of great importance for the security of c.

VII. GOING MULTI-PARTY

To be able use the transformation in a multi-party scenario, we propose to use secure computations based on Shamirshared values as developed by Ben-Or et al. [START_REF] Ben-Or | Completeness theorems for non-cryptographic fault-tolerant distributed computation[END_REF]. In short, this technique allow us to share values among the p parties in a way that the knowledge of less than k < p/2 of the p shares does not reveal any information about the secret value. Yet these shares can be used to make computations (additions, multiplications etc.) on the secret values. After finishing the computations, we can put the shares of the result together and reconstruct it -whereas all input and intermediate values remain secret.

As in our scenario the input values (M {1,2} , b {1,2} , c) are distributed among the parties, we have to assemble the data and jointly choose the random values (P , Q, S, A, r). The assembly of the data depends on the initial partitioning in the problem. However -as long as each value is clearly "owned" by one party (which is usually the case) -the assembly is not a problem: Each party shares its values and the matrices and vectors are filled with the shares.

When choosing the random matrices and vectors, we have to ensure the randomness and the correct form of the result, i.e. guarantee that Q is a monomial matrix, A a permutation matrix and so on. Moreover we need to ensure that no party knows the resulting values, because they would allow the parties to undo the transformation and access the secret data.

The basic idea is that each party randomly chooses its random values (i.e. party i chooses matrices P i , Q i , S i , A i and a vector r i ) and we combine them to have jointly and fairly chosen random values. Ideally, even the knowledge of p -1 of these p input values should not leak any information about the result. To achieve this, we rely on secure computations again. After each party has chosen its input, we calculate

P = p-1 i=0 P i , Q = p-1 i=0 Q i , S = p-1 i=0 S i , A = p-1 i=0 A i and r = p-1 i=0 r i .
As these computations take place in a finite field, the result is random, and even knowing p -1 of the p input values does not reveal anything about the output. This leads to following protocol: 1) Each party i chooses a random invertible matrix P i , a random positive monomial matrix Q i , a random positive diagonal Matrix S i , a random permutation Matrix A i and a random positive vector r i . 2) Each party i securely shares its parts of M 1 , M 2 , b 1 , b 2 and c as well as P i , Q i , S i , A i and r i . 3) Secure computation of P , Q, S, A and r as described above. Assembly of M 1 , M 2 , b 1 , b 2 and c. Transformation:

Calculation of M = P *   M 1 Q 0 M 2 Q A -S 

A. Optimization: Fast Matrix Multiplication

Matrix multiplication is the key operation of our transformation. When using secure computations on Shamir-shared values, multiplications are expensive due to the degree-reducing step on the shares. We want to reduce this overhead as much as possible.

A naive implementation (see Algorithm VII.1) of the multiplication of a shared (m × k) matrix A and a (k × n) matrix B will result in a time complexity of O(mkn), a round complexity (i.e. the number of points in time when communication is necessary) of O(mkn) and a communications complexity (i.e. the number of messages exchanged) of O(mknp 2 ) as each degree-reducing step requires one round and O(p 2 ) messages [START_REF] Ben-Or | Completeness theorems for non-cryptographic fault-tolerant distributed computation[END_REF]. This can be reduced to one round and O(mnp 2 ) messages by postponing the degree-reducing step as much as possible.

If A x,y and B y,z are shares of a polynomial with degree k, then temp is of degree 2k. Yet, as adding two shares does not increase the degree of the polynomial, we can calculate C x,z = k-1 y=0 A x,y * B y,z without exceeding a degree of 2k. Thus we can postpone the degree-reducing step until the end of all calculations and execute it in parallel on the whole matrix to reduce round complexity to 1 and communication complexity to O(mnp 2 ) (see Algorithm VII.2 6 ). for z = 0 to n -1 do end for 9: end for 10: execute degree-reducing step on all m * n entries of C in parallel VIII. USE CASE: SUPPLY CHAIN OPTIMIZATION Supply Chain Optimization problems can be expressed as linear programs [START_REF] Stadler | Supply chain management and advanced planning[END_REF]. Each variable corresponds either to the number of units to produce at a certain factory at a certain time, to the inventory (per factory and time), or to transported units (between two factories). The constraints are flow constraints (e.g. "all produced units are either shipped to the next stage in the chain or remain in the inventory"), capacity constraints (e.g. "factory X can produce at most n units of product Y in period i") and demand constraints (e.g. "in period i, k units of product Z are sold to customers"). The cost function is composed of production costs (per produced unit), shipping costs (price of shipping one unit from factory X to factory Y ) and inventory costs (the price of stocking products or intermediate components at a factory). This means that variable is "owned" either by one company (for production an inventory variables) or by two companies (for transport variables). Hence we have a simple structure which allows us to easily set up the distributed linear program as discussed above. Furthermore, most of the constraints (the flow and demand constraints) are expressed as equations, which is convenient for our transformation.

A. Experimental Results

To analyze the practical performance of our approach, we implemented a prototype. The test data is a sample supply chain structure with five companies. To scale the size of the problem, this optimization was done for different numbers of planning periods which increases the number of variables and constraints. The measured timings include sharing of the values, performing the transformation, reconstructing the transformed LP, solving it using a standard LP solver and backtransforming the solution using secure computation again (for more details see the technical report [START_REF] Dreier | Practical Secure and Efficient Multiparty Linear Programming Based on Problem Transformation[END_REF]).

The results (see Figure VIII-A) are promising: The standard test case (six periods: 180 constraints, 282 variables) is solved in approximately 25 minutes, and even the bigger test cases are completed within two hours. This is sufficiently fast for practical usage and much quicker than the distributed Simplex algorithms.

Based on the bounds established above, we calculated the numerical values for the a-posteriori probabilities of a right guess in the six-period case for a uniform input distribution: As explained above, there is no easily calculable bound on the leakage of the constraints, but the chances of guessing P only based on A are less than 2.34 • 10 -746370 .

IX. CONCLUSION

In this paper we proposed a disguising transformation for linear programs. We proved correctness and analyzed security in the framework of leakage quantification by establishing bounds on the leakage of the secret data. Subsequently we developed a secure multi-party protocol to make the transformation usable in a distributed case.

Based on this protocol we implemented a prototype for secure supply chain optimization. We conducted experiments which show significant performance gains compared to entirely distributed secure simplex algorithms.

Moreover, the basic ideas of the transformation can easily be extended and generalized to fit other problems such as systems of linear equations. This highlights the power of disguising approaches.

A. Future Work

We can imagine several extensions to the transformation to enhance security, for example splitting variables or creating fake variables in order to protect c. These measures have to be checked for their security as well as their impact on performance.

Furthermore the numerical stability of the transformation should be examined. During our experiments a standard LP-Solver was always able to find the correct solution, however rounding errors and performance deterioration could be observed from time to time. This also is of interest as numerical effects could lead to attacks.

The security analysis could also be refined in some parts. In particular the interconnections between the different parts of the transformation should be examined more precisely as this is one of the most probable entrance point for attacks. Additionally it would be desirable to refine the bounds on the leakage of combined channels (i.e. to show that for example the combination of a multiplication and a permutation is strictly better than the best of them alone) and give -if possible -a closed-form bound on the leakage of a matrix multiplication.

As our security analysis currently supposes a Honest-but-Curious-Scenario, it could also be interesting to analyze the possible influence of malicious participants with byzantine behavior (i.e. participants that give wrong input data or do not follow the protocol) on the results. In the case of distributed supply chain management, it would be particularly interesting to examine if parties are able to obtain a financial advantage by manipulating the outcome using specially prepared input data, and if this would be noticeable to the other parties.

Fig. 1 .

 1 Fig. 1. Privacy-Preserving Linear Programming based on Problem Transformation

b z s ≥ 0 6 )

 6 = c T * Q 0 . . . 0 . 4) c T s , M and b are reconstructed from the shares and passed into the cloud. 5) The Cloud solves the Linear Program min c T s z s s.t. M z s = Secure sharing of the solutions z 7) Secure distributed computing of x = Q(z -r) 8) Output of values of x to their respective owner.

Algorithm VII. 1 2 : 3 :← 0 4 : 2

 12342 Matrix Multiplication C = A * B, naive 1: for x = 0 to m -1 do for z = 0 to n -1 do sum for y = 0 to k -1 do 5: temp ← A x,y * B y,z Matrix Multiplication C = A * B, optimized communication 1: for x = 0 to m -1 do 2:

4 :

 4 for y = 0 to k -1 do 5: sum ← sum + A x,y * B y,z

Fig. 3 .

 3 Fig. 3. Performance in Distributed Supply Chain Management Case

•

  cost function c: P R posteriori ≤ 3.76 • 10 -220 • values x: P R posteriori ≤ 2.30 • 10 -1409

A monomial matrix contains exactly one non-zero entry per row and column.

A diagonal matrix where the entries on the main diagonal A i,i are strictly positive.

A permutation matrix is a monomial matrix where the non-zero entries are "1".

A slack-variable is used to express inequality constraints as equality constraints. The idea is to introduce an additional variable for each constraint which takes up the "remainder" or "slack". For example, instead of 3x 1 ≤ 10 we can write 3x 1 + s 1 = 10 for s 1 ≥ 0.

Braun et al. [5] call this value the Multiplicative Leakage. The interesting point is the maximum of this value over all X, which we call Multiplicative Leakage.

This algorithm is somewhat similar to the inner product algorithm of Catrina and de Hoogh[START_REF] Catrina | Secure multiparty linear programming using fixed-point arithmetic[END_REF], but generalized to a complete matrix multiplication.