M Ciesielski

S Askar
email: saskar@ecs.umass.edu

D Gomez-Prado

J Guillot
email: jguillot@univ-ubs.fr

E Boutillon
email: emmanuel.boutillon@univ-ubs.fr

Data-Flow Transformations using Taylor Expansion Diagrams

An original technique to transform functional representation of the design into a structural representation in form of a data flow graph (DFG) is described. A canonical, word-level data structure, Taylor Expansion Diagram (TED), is used as a vehicle to effect this transformation. The problem is formulated as that of applying a sequence of decomposition cuts to a TED that transforms it into a DFG optimized for a particular objective. A systematic approach to arrive at such a decomposition is described. Experimental results show that such constructed DFG provides a better starting point for architectural synthesis than those extracted directly from HDL specifications.

I. INTRODUCTION

Although considerable progress has been made in behavioral and high-level synthesis over the past two decades, new solutions are needed for advanced algorithm-oriented designs, such as multi-media and signal processing applications. Traditional high-level synthesis tools are effective at capturing the HDL specification by extracting a data flow graph (DFG) and mapping it directly into an architecture using structural optimization techniques, such as scheduling, resource allocation and binding [START_REF] Demicheli | Synthesis and Optimization of Digital Circuits[END_REF]. However, these tools are less effective at the higher levels of abstraction. Most of these systems rely on a fixed data flow graph (DFG), extracted from the initial HDL specification, and provide limited means for modifying the initial data flow structure. In an attempt to explore other solutions the user often needs to rewrite the original specification, from which another DFG is derived and synthesized. This approach seriously limits the scope of the ensuing architectural synthesis and the quality of final hardware implementation.

Several attempts have been made to provide optimizing transformations in high-level synthesis [START_REF] Guerra | High level Synthesis for Reconfigurable Datapath Structures[END_REF], [START_REF] Chaiyakul | High-level transformations for minimizing syntactic variances[END_REF], [START_REF] Chandrakasan | Optimizing Power Using Transformations[END_REF], [START_REF] Potkonjak | Synthesis for testability using transformations[END_REF], [START_REF] Srivastava | Optimum and heuristic transformation techniques for simultaneous optimization of latency and throughtput[END_REF], [START_REF] Püschel | SPIRAL: Code Generation for DSP Transforms[END_REF], [START_REF] Gupta | Using Global Code Motion to Improve the Quality of Results in High Level Synthesis[END_REF]. Behavioral transformations have been also used in other areas, such as optimizing compilers [START_REF] Ullman | Computational Aspects of VLSI[END_REF] and logic synthesis [START_REF] Demicheli | Synthesis and Optimization of Digital Circuits[END_REF]. With the exception for a few specialized systems for DSP code generation, such as SPIRAL [START_REF] Püschel | SPIRAL: Code Generation for DSP Transforms[END_REF], these methods rely on simple manipulations of algebraic expressions based on term rewriting and algebraic properties of associativity, commutativity, and distributivity. Several high-level synthesis systems, such as Cyber [START_REF] Wakabayashi | Cyber: High Level Synthesis System from Software into ASIC[END_REF], Spark [START_REF] Gupta | SPARK: A Parallelizing Approach to the High-Level Synthesis of Digital Circuits[END_REF] and others, use a host of methods for code optimization (such as kernelbased algebraic factorization, branch balancing, speculative code motion methods, dead code elimination, etc.) but without relying on any canonical representation that would guarantee local optimality of the transformations. To the best of our knowledge no systematic method for DFG modification for the purpose of behavioral optimization and synthesis has been presented to date.

In contrast, this paper presents a systematic method to perform behavioral transformations, based on a canonical representation of the computation, called Taylor Expansion Diagram (TED) [START_REF] Ciesielski | Taylor Expansion Diagrams: A Canonical Representation for Verification of Data Flow Designs[END_REF]. Being a canonical representation, TED can capture an entire class of structural solutions, rather than a single DFG. Through a unique decomposition procedure proposed in this paper, this functional representation is converted into a structural representations (DFG), optimized for a particular design objective. This approach provides means for fast design space exploration directly from behavioral specifications.

A. Canonical TED Representation

Taylor Expansion Diagram [START_REF] Ciesielski | Taylor Expansion Diagrams: A Canonical Representation for Verification of Data Flow Designs[END_REF] is a canonical, word-level data structure that offers an efficient way to represent computation in a compact, factored form. It is particularly suitable for algorithm-oriented applications, such as signal and image processing, with computations modeled as polynomial expressions.

An algebraic, multi-variate expression, f (x, y, ...), can be represented using Taylor series expansion, w.r.t. variable x as follows:

f (x, y, . . .) = f (x = 0) + xf ′ (x = 0) + 1 2 x 2 f "(x = 0) + . . .
where f ′ (x), f ′′ (x), etc, are the successive derivatives of f w.r.t. x. The terms of the decomposition are then decomposed with respect to the remaining variables (y, .., etc.), one variable at a time. The resulting decomposition is stored as a directed acyclic graph whose nodes represent the terms of the expansion. The number of children at each node depends on the order of the polynomial expression (w.r.t its decomposing variable) rooted at that node. The resulting graph is reduced, normalized and canonical for a fixed ordering of variables. The expression represented by the graph is computed as a sum of the expressions of all the paths of the graph, from root to terminal 1.

Figure 1(a) shows one-level decomposition of function f (x, y, . . .) at variable x. The nodes f (x = 0, y, . . .), f ′ (x = 0, y, . . .), etc, represent subsequent derivative functions that depend only on the remaining variables: y, etc. For the ordering of variables (A, B, C), the reduced and normalized TED is constructed as follows. The decomposition is performed first with respect to variable A, producing the following terms:

F (A = 0) = 2BC, F ′ (A = 0) = B + 2C, and 1 2 • F ′′ (A = 0) = 1.
Next the expansion is applied to the resulting non-trivial terms F (0), F ′ (0) with respect to variable B, and subsequently with respect to variable C. Note the multiplicative weights assigned to the edges (default weight is 1).

In summary, TED is a multiplicative diagram which maps word-level (integer) inputs into integer, word-level (integer) outputs. The function encoded in the TED is computed as a sum of the expressions of all the paths from root to terminal node 1. (For this reason, TED can be viewed as a graph without edges leading to terminal 0, and containing only terminal node 1.) The expression of each path is computed as a product of the expressions of all the edges in the path; the expression of each edge is in turn a product of the variable in its respective power, weighted by the numerical label (coefficient) assigned to the edge. For example, the expression for the function encoded in the TED in Figure 2(a) is computed as a sum of two paths from root to terminal node 1:

F = A•B+A•B 0 •C = A•B+A•C, or equivalently F = A•(B+C).
The latter form illustrates the fact that TED encodes useful factorization of the expression. This can be seen by noting that subexpression A, corresponding to the linear edge spanning out of variable A, is common to the two paths, A • B and A • C, and can be factored out from the expression. This factorization is manifested in the graph by the presence of a subexpression (B + C), rooted at node B, which can be extracted from the graph. This is an important feature of the TED representation, explored in factorization and common subexpression extraction [START_REF] Guillot | Efficient Factorization of DSP Transforms using Taylor Expansion Diagrams[END_REF].

In this context TED can serve as a canonical view of the computation, regardless of the specific data flow or architectural implementation. Several structural, data flow solutions can be derived from such a representation using the transformation process described here, guided by the desired objective (area, latency or power), resulting in an efficient hardware implementations. and requires two multipliers and one adder, while the one in Figure 2(c) minimizes the number of operators at a cost of increased latency (L = 3). The two solutions correspond to different scheduling of the same DFG. However, another, better solution can be obtained by transforming the original specification

B. Motivation

Consider a simple computation,

F = A • B + A • C,
F = A • B + A • C into F = A • (B + C).
The modified expression corresponds to a different DFG, which uses only one adder and one multiplier, and which can be scheduled in two control steps (L = 2), as shown in Figure 2(d). The existing commercial tools do not perform such transformations; they only derive a single, fixed data flow graph (DFG) for the given computation and have no capability to transform it into another DFG (better from a given objective point of view) prior to the architectural synthesis. In order to obtain an optimum solution (or to explore other solutions), the user must rewrite the original specification, from which another fixed DFG is derived and subjected to the architectural optimization.

To address this problem we propose a systematic way to transform the initial design specification into a data flow graph, optimized for a given objective. Figure 2 shows the initial HDL specification represented as a canonical TED along with all the solutions that can be obtained from this canonical representation. In particular, Figure 2(d) shows the solution with minimum latency subject to the imposed resource constraint.

II. BEHAVIORAL TRANSFORMATIONS

A. System Flow

The basic idea behind the proposed system is to transform the functional TED representation of the design into a structural DFG representation. A partially scheduled DFG is obtained from TED by performing successive decomposition of the TED by means of cuts. The cut-based decomposition is guided in such a way as to optimize the DFG for a given objective. Input to the system is an algorithmic description of the design in C or behavioral HDL. The system flow, shown in Figure 3, is composed of the following steps. High-level synthesis system, GAUT [START_REF] De | GAUT, Architectural Synthesis Tool[END_REF], is used for front-end compilation and architectural synthesis, while the selection of best DFG is performed by cut-based TED decomposition.

1) Compilation: Compile the design to generate functional representation (TED). Minimize TED using variable ordering [START_REF] Gomez-Prado | Variable Ordering for Taylor Expansion Diagrams[END_REF]. This phase is intended to minimize the number of ADD/MULT operations in the resulting DFG. 2) TED Factorization and Extraction: Generate hierarchical TED representation by performing TED factorization and common subexpression extraction using techniques described in [START_REF] Guillot | Efficient Factorization of DSP Transforms using Taylor Expansion Diagrams[END_REF]. Hierarchical TED is a graph, whose nodes represent TEDs of local functions. 3) Transformation of TED into DFG: Perform decomposition of local TEDs using sequence of cuts to create a partially scheduled DFG. The decomposition process is guided by the desired objective, such as latency, operators area, etc.

This paper concentrates on step 3 of the flow, using latency as a primary objective. Figure 4 shows a TED for function P = (G + H) + F (I + J). (To simplify the explanation, TED is represented with a single terminal 1, with all the edges leading to terminal node 0 removed.) Cut A3 in the figure partitions function P disjunctively into (G + H) and F (I + J). This is shown as the top level decomposition in Figure 5. Cut M 1 is an example of a multiplicative cut. It breakes the right TED subgraph, Q = F (I + J), generated by applying cut A3, into two expressions: F and (I+J) to be composed conjunctively (refer to Figure 5). Each time an additive or multiplicative cut is applied to a TED, a hardware operator (ADD or MULT) is introduced in the DFG to perform the required operation on the two subexpressions. This way, a functional TED representation (composed of algebraic operations) is eventually transformed into a structural representation (composed of hardware operators), a data flow graph (DFG).

C. Decomposition Algorithm

The additive (multiplicative) cut is called admissible if it partitions the TED expression into exactly two subexpressions in an additive (multiplicative) way, and it does not increase the number of operations determined by the original TED representation. Otherwise the cut is called inadmissible. Note that admissibility is a dynamic property; a cut may be inadmissible for a T ED(k) in step k, and become admissible for one of the subgraphs of T ED(j), at some step j > k. For example, cut M 1 in Figure 4 is admissible only after the cut A3 has been applied. Subsequent application of M 1 makes cut A2 admissible, etc. Finding an admissible cut sequence is one of the main tasks of this work. A cut-based decomposition is illustrated in Figure 5 for an admissible sequence (A3, M 1, A2, A1). Lemma: Each admissible cut sequence generates a unique DFG.

Proof: Each sequence, being admissible, produces a legal DFG. Furthermore only one such DFG can be obtained for a given sequence since different DFGs mean different partial order of DFG nodes (operators), hence different cut sequence.

The resulting DFGs may differ in terms of critical path, balance structure, etc., which translates into physical characteristics of the designs (latency, power dissipation due to spurious computation, etc.), synthesized from these DFGs. An ordering of cuts is sought that will optimize the DFG structure for a given design objective (tree balance, longest path, etc).

We use a Branch & Bound algorithm to search for admissible cut sequences that optimize a particular design objective. Here we concentrate on minimizing latency, i.e., the critical path of the expected DFG. The algorithm uses cut admissibility for branching, and latency of the current solution as the bounding condition.

III. RESULTS

The TED-to-DFG transformation procedure has been implemented as part of an experimental tool TEDify.

[16], available on the web. The generated DFGs have been compared to those obtained by the GAUT high-level synthesis system [START_REF] De | GAUT, Architectural Synthesis Tool[END_REF] for a number of DSP designs.

Figure 7 shows the DFG extracted by GAUT from the original VHDL or C description of an FIR16 filter. This DFG is then used in subsequent architectural optimization. Figure 8 shows the DFG obtained by our TEDify tool: a TED was extracted from the same specification and was subjected to a cut-based TED decomposition to produce a DFG with minimum latency.

The reduction in the number of multiplier operations, from 16 required by the original specification extracted by GAUT to 8 produced by TEDify, has been accomplished by the first two steps of the flow (1: TED variable ordering; and 2: TED factorization and expression extraction), while the number of adders remained the same. The latency of a DFG produced in Table I shows the pre-synthesis results; it compares the DFGs (in terms of the number of operators and control steps) produced by GAUT and TEDify for a number of DSP designs (written in C or behavioral VHDL). In most cases the latency was reduced, sometimes drastically. GAUT AND DFG GENERATED BY TEDIFY Figure 9 shows the Gantt chart for the FIR-16 filter, synthesized by GAUT using the original DFG (extracted by GAUT), with the number of resources limited to one ADD and one MULT. (see row 1 of Table II). Figure 10 shows the Gantt chart for the same design under the same resource constraint, synthesized by GAUT using the DFG generated by TEDify. We can see the difference in latency (210 ns vs 330 ns) and in the number of registers needed by the two solutions. We are currently using TED ordering which minimizes the number of nodes of the TED, but other objectives that will more directly address the number of operators (especially the multipliers) should be investigated. The correlation between the number of multiplicative edges and multipliers would suggest that ordering which minimizes the number of those edges, rather than the number of TED nodes, would be more appropriate for resource-limited applications. This objective should also drive the factorization and extraction algorithms, which are an important part of the TED decomposition scheme, prior to cut-based transformation onto DFG. Those and other practical issues are currently under investigation.

Furthermore, the area cost should be properly qualified by distinguishing multipliers with arbitrary operands from those that use constants, especially if the constants are powers of 2, in which case the multiplier can be replaced by a less expensive shifter.

Finally, there are natural limitations inherent to the TED representation. Currently TEDs can handle designs with adders, subtractors, multipliers, shifters and Boolean logic, but cannot efficiently handle other operators. For this reason the use of TEDs is currently limited to untimed algorithmic descriptions of DSP designs.

2 Fig. 1 .Figure 1

 211 Fig. 1. TED: a) Decomposition principle; b) TED example for F = A 2 + AB + 2AC + 2BC Figure 1(b) shows TED for function F = A 2 + AB + 2AC + 2BC.Additive edges, corresponding to F (0) are represented in the graph as dotted lines. Linear (first order) edges, associated with F ′ (0) are shown as solid lines. The second order edges, associated with F ′′ (0) are shown and double lines. For the ordering of variables (A, B, C), the reduced and normalized TED is constructed as follows. The decomposition is performed first with respect to variable A, producing the following terms:F (A = 0) = 2BC, F ′ (A = 0) = B + 2C, and 1 2 • F ′′ (A = 0) = 1.Next the expansion is applied to the resulting non-trivial terms F (0), F ′ (0) with respect to variable B, and subsequently with respect to variable C. Note the multiplicative weights assigned to the edges (default weight is 1).

  where variables A, B, C are word level signals. The TED representation for this computation is shown in Figure 2(a).

Figure 2 (

 2 b) and (c) show two architectural solutions that can be obtained using state-of-the art commercial synthesis tools. The solution shown in Figure 2(b) minimizes latency (L = 2)

Fig. 2 .

 2 Fig. 2. Behavioral transformations: a) Canonical TED representation; b,c,d) Functionally-equivalent DFGs.

Fig. 3 .

 3 Fig. 3. System flow: behavioral transformation and DFG generation provided on top of architectural synthesis

Fig. 4 .

 4 Fig. 4. Cuts in TED for P = G + H + F (I + J).

Fig. 5 .

 5 Fig. 5. Decomposition based on cut sequence (A3, M 1, A2, A1) on TED for P in Figure 4.

Figure 6 Fig. 6 .

 66 Figure6shows two DFGs derived from two different cut sequences, (A3, A1, M 1, A2) and (A1, A3, M 1, A2), applied to the TED in Figure4.The following lemma is an obvious consequence of the cutbased procedure described above.

Fig. 7 .

 7 Fig. 7. DFG for FIR16 produced by GAUT

Fig. 9 .

 9 Fig. 9. Gantt chart for FIR16 design synthesized with the original DFG extracted by GAUT

Table

  II shows the results for FIR-16 and IIR filter designs synthesized by GAUT under resource constraints for the original DFG, extracted by GAUT and for the DFG computed by TEDify. The resource limits imposed on the number of ADD and MULT operators are shown in columns 2 and 3 of the table. The next two columns show the latency (in ns) of the final architecture for the two DFGs (produced by GAUT

	x5	x10	x6	x9	x7	x8	xn	x15	x4	x11	x1	x14	x13	x2	x12	x3
	-2097	+ 5	-1674	+ 6	9883	+ 7	98	+ 8	950	+ 4	-39	+ 3	+ 1	-327	+ 2	439
		* 5	* 6			* 7	* 8			* 4	* 3		* 1		* 2
			+ 11			+ 12				+ 10			+ 9
						+ 14				+ 13
								+ 15
								+ yn
	Fig. 8. DFG for FIR16 produced by TEDify
						DFG: GAUT			DFG: TEDify
			Design		# Operations	Latency	# Operations	Latency
					ADD	MULT	c-steps	ADD		MULT	c-steps
			IIR		4	4		5		4		0		3
			FIR16		15	16		16		15		8		5
			Prodmat		48	64		4		48		64		3
			Ellip		20	0		11		36		0		7
			DCT		48	64		4		48		64		3

TABLE I PRE

 I -SYNTHESIS RESULTS: COMPARISON OF DFG EXTRACTED FROM THE INITIAL SPECIFICATION BY GAUT AND DFG PRODUCED BY TEDIFYand TEDify, respectively), synthesized by GAUT. The delay of ADD was set to 10 ns and MULT to 20 ns, with the clock cycle set to 10 ns (these parameters are user-defined and can be altered).

	Benchmark	# of Operators	Latency(ns)
		ADD	MULT	GAUT	TEDify
	FIR16	1	1	330	210
		1	2	180	150
		1	4	170	150
		1	8	170	150
		2	1	330	210
		2	2	180	130
		2	4	170	90
		2	8	170	90
		4	1	330	210
		4	2	180	130
		4	4	170	90
		4	8	170	70
	IIR	1	1	90	40
		1	2	90	40
		2	1	60	30
		2	2	60	30

V. ACKNOWLEDGMENT

This work has been supported by a grant from the National Science Foundation under award No. CCR-0204146.


O Franco

X Beudaert
email: xbeudaert@ideko.es

K Erkorkmaz

A Barrios

I Mancisidor

J Munoa

Control strategies for structural chatter suppression using a rack and pinion drive

Keywords: Machine tool, control, chatter, drive, active damping 1

Chatter vibrations are a major limitation for ma chining operations and lead to productivity reduction, low tool life and poor surface finish. Active damping with the machine tool's own drives can be used to increase the stability limit related to structural modes of the machine. Various feedback strategies can be used measuring the vibrations with an external accelerometer and feeding back this signal to the position, velocity or current setpoint of the drive. The objective of this article is to build a model based on experimental data to predict the effect of the standard control parameters and additional acceleration feedback controller on the Frequency Response Function of the machine. Experimental tests are carried out on an electronically preloaded rack and pinion drive. The open loop Frequency Response Functions of the dynamical system are measured and the closed loop responses are calculated using Linear Fractional Transformations representing the control structure of the Multi Input Multi Output system.

Introduction

In heavy duty milling operations, machine tool users are looking for high material removal rates which can be limited by chatter vibrations [START_REF] Munoa | Interaction between multiple modes in milling processes[END_REF]. For structural chatter vibrations coming from the low frequency modes of the complete machine it has been shown that the use of active techniques is justified because they can cope with the variations of the machine dynamics according to the varying cutting position [START_REF] Munoa | Chatter suppression techniques in metal cutting[END_REF]. For heavy duty machining operations, powerful actuators have been designed using hydraulic [START_REF] Brecher | Electrohydraulic Active Damping System[END_REF] or electromagnetic technologies [START_REF] Munoa | Chatter suppression in ram type travelling column milling machines using a biaxial inertial actuator[END_REF]. The use of an additional actuator close to the cutting point increases significantly the cost and requires a modification of the machine design. Instead of installing a new actuator, it is possible to use the machine's own drives to actively damp the low frequency modes.

Active damping techniques using machine feed drives have been applied mostly to reduce inertial vibrations coming from the machine movements. A comprehensive review is presented by Altintas et al. in [START_REF] Altintas | Machine tool feed drives[END_REF]. Hosseinabadi and Altintas [START_REF] Hosseinabadi | Modeling and active damping of structural vibrations in machine tools[END_REF] modelled a ball-screw drive to simulate the effect of the additional active damping loop in time domain and to tune a sliding mode controller. Chatter suppression using machine drives has been simulated first by Chen and Tlusty [START_REF] Chen | Effect of low-friction guideways and lead-screw flexibility on dynamics of high-speed machines[END_REF], and then it has first been applied by Munoa et al. [START_REF] Munoa | Active suppression of structural chatter vibrations using machine drives and accelerometers[END_REF]. A significant increase in the chatter stability limit is demonstrated experimentally on a large milling machine equipped with a rack and pinion drive. Later, the limitations of this chatter suppression technique has been analysed to emphasize the constraints related to the non-collocated control, to the technical implementation, and to the drive bandwidth [START_REF] Beudaert | Limiting factors for the active suppression of structural chatter vibrations using machine's drives[END_REF]. The importance of the standard cascade control parameters has been remarked however, up to now, their interaction with additional active damping controller has not been studied in details.

For medium and large size heavy duty turning / milling machines, the rack and pinion, electrically preloaded, drive solution is commonly used to handle the long travel and high loads [START_REF] Uriarte | Machine tools for large parts[END_REF]. As opposed to the linear and ball-screw drives which concentrated most of the attention of the research community, the rack and pinion drives have seldom been studied in the literature. Zirn [START_REF] Zirn | Machine Tool Analysis-Modelling, Simulation and Control of Machine Tool Manipulators[END_REF] analysed the rack and pinion drives from the control point of view. Engelberth et al. [START_REF] Engelberth | Properties of electrically preloaded rack-and-pinion drives[END_REF] analysed the effect of the electronic preload on the static stiffness, bandwidth and backlash.

The objective of this paper is to predict the influence of the control on the closed loop frequency response functions using measured Frequency Response Functions on a rack and pinion drive. Knowing the position and velocity loops control structure, the effect of the gains of the P-PI controllers can be simulated as well as the introduction of the active damping loops.

The rest of the paper is organized as follows. The model of a preloaded rack and pinion drive is presented in Section 2. Then, Section 3 presents the multi input multi output frequency domain analysis of the system. Section 4 contains the experimental validation and the active damping demonstration. Finally, this paper is concluded in Section 5.

2.

Modelling of rack and pinion drives Figure 1 shows the control scheme of a rack and pinion drive with a Master-Slave control structure. The position loop is closed only for the master motor based on the linear scale measurement to compute the velocity setpoint set  . Each motor is controlled by a PI velocity controller according to the measured velocity on the motor encoders ( 1 and 2   ) and generate a torque ( 1  and 2  ). The current loops are neglected as they have no influence at low frequency. To apply the electronic preload, the torque preload t M is applied using a third PI controller which generates a non-uniform distribution of the torque reference ref  on both motors. For two identical motors, the weighting factors TM K and TS K have to be chosen as equal (e.g. 50%).

When the drive is stopped ( 0 ref  

), both motors generate a torque in opposite directions which ensures that the pinions are in contact with different sides of the rack. This value of the torque preload t M is commonly specified as a percentage of the nominal torque of the motor, typically in the range of 10 to 30% [START_REF] Zirn | Machine Tool Analysis-Modelling, Simulation and Control of Machine Tool Manipulators[END_REF]. Additionally to the closed loop controls, the velocity and acceleration feedforward actions can be used to decrease the tracking errors. The machine studied in this paper is using a velocity feedforward but is has no effect on the machine tool's feedback dynamics. Hence, in the forthcoming analysis, its contribution can be neglected.

On top of this standard control architecture, three additional commands are summed respectively to the position, velocity and torque setpoints ( d x , d  , d  ) to apply the active damping based on the measurement of the ram tip acceleration. The objective of the active damping loops is to minimize the ram tip compliance 2 / x F c   .

Frequency Domain Analysis

The objective is to analyse the effect of the control of the response of a rack and pinion drive purely in frequency domain using experimental data. Hence, the tedious task of model fitting is not required and all the complexity of the machine dynamics is implicitly contained in the measured frequency response function data.

Multi input multi output model

The general control scheme of a rack and pinion drive shown in Figure 1 is modelled using the Linear Fractional Transformation technique (LFT) (see Figure 2). The Linear Fractional Transformation is a convenient and powerful formulation for MIMO systems, it allows decomposing the whole system into open and closed loops signals to work independently with each part. The block which defines the open loop behaviour of the rack and pinion drive is defined as matrix P which is here expended by set  and t M to satisfy the mathematical requirements of a Linear Fractional Transformation. The speed controller is defined as matrix K composed of 2x4 Frequency Response Functions (FRFs). With the formalism of the LFT, the closed velocity loop frequency response H is obtained in Eq. 1 [START_REF] Skogestad | Multivariable Feedback Control: Analysis and Design[END_REF]. (1)

Discretization of the speed controller

Using the speed controller structure highlighted in Figure 1, the eight transfer functions between the inputs variables defined in the matrix U and the torque commands (matrix Z) can be obtained. The relationships for the master drive's torque command are described in the Eq. (2-5), similar expression are derived for the slave's torque command.

             2 S T M T M T 2 S T S T S T 2 M 2 T S T 1 T 2 2 S T 1 1 TS TM T S TS s Ti M Ti Tn s S Ti M Ti s Ti Tn Ti K S Kp K s Ti M Ti Tn s S Ti M Ti s Ti Tn Ti K M Kp K s S Ti Ti Tn s s Ti Tn s Ti s Ti K S Kp K s Ti Ti s M Ti Kp D M               (2)        2 S T M T M TS T 2 S T S T S TM T 2 M T 1 1 2 M T T s Ti M Ti Tn s S Ti M Ti s Ti Tn Ti K S Kp K s Ti M Ti Tn s S Ti M Ti s Ti Tn Ti K M Kp K s S Ti Ti Tn s Ti s Tn Ti Kp K D S M            (3)               2 S T M T M TS T 2 S T S T S TM T 2 M T 2 S T S T 1 T 2 S T 1 3 TS s Ti M Ti Tn s S Ti M Ti s Ti Tn Ti K S Kp K s Ti M Ti Tn s S Ti M Ti s Ti Tn Ti K M Kp K s S Ti Ti Tn s s Ti Tn s Ti s Tn K S Kp Kp s Ti Tn s M Ti M Kp D               (4)           2 S T M T M TS T 2 S T S T S TM T 2 M T 1 1 1 TS T 4 s Ti M Ti Tn s S Ti M Ti s Ti Tn Ti K S Kp K s Ti M Ti Tn s S Ti M Ti s Ti Tn Ti K M Kp K s S Ti Ti Tn s s S Ti s M Ti s T Tn K S Kp M Kp K D            
(5)

Position and active control loops

As discussed above, there are different control loops where the active control can influence. addition of the torque command d  is implemented by modifying the velocity loop LFT such that the controller K also takes the ram tip acceleration 2 x  as input. Hence, the action of the damping controller is added to the action of the speed controller. The position loop, the velocity and position command damping are modelled by another Linear Fractional Transformation considering the closed velocity loop H obtained previously. The controller of the position and velocity command damping loops are dp K and dv K , respectively (see Figure 2b).

Experimental validation

Description of the experiments

The experiments are realized on a large vertical turning center of SORALUCE equipped with a SIEMENS 840D Solution Line CNC, see Figure 3a. Two fast analogue inputs are connected to the Profinet bus of the CNC to measure the ram tip acceleration  is used to apply a chirp excitation to the drive and to obtain the P12 and P22 matrices of FRFs. Then, the position and velocity loops gains are reduced to zero and a chirp excitation is applied by a shaker to obtain the FRFs P 1,1 , P 2,1 , P 5,1 , P 6,1 depending on the force input c F . For rack and pinion drives, the electronic preload plays an important role as seen in [START_REF] Engelberth | Properties of electrically preloaded rack-and-pinion drives[END_REF]. All the tests are realized with 10% preload not to change the open loop behaviour of the system.

Validation of the predictions

To validate the predictions, the measured and predicted ram tip FRFs are shown in Figure 3b. The Linear Fractional Transformation could predict the tendency to increase the ram tip compliance when the proportional gain of the velocity loop Kp is increased. However, the accuracy of the prediction is not satisfying to rely on this model to tune the control parameters. Further analyses are required to handle the non-linear open loop response of the system and to quantify the impact of the non-linearity on the FRF predictions.

Active damping with position, velocity or acceleration additional loop

Figure 4a presents the Frequency Response Function at the ram tip with and without active control obtained with an impact hammer. It can be seen that the three options using position, velocity or acceleration additional commands can provide an important reduction of the compliance amplitude. A more detailed performance evaluation is required to determine which option should be preferred.

Finally, the vibration amplitude measured during a turning cutting test is shown in Figure 4b. It can be seen that the vibration amplitude is decreased when the active control is applied. However, chatter vibration still persists and further improvements are required to be able to significantly increase the chatter stability limit.

Conclusion

This paper has presented a method to analyse the behaviour of a rack and pinion drive in the frequency domain. Three additional control loops can be added to the standard control structure to actively damp the low frequency modes of the machine structure. The Linear Fractional Transformations presented here allow predicting the effect of the control on the ram tip compliance. However, further investigations are needed to get accurate predictions. Indeed, the complex behaviour of the preloaded rack and pinion drive should be evaluated more accurately to build a reliable model of the influence of the control on the dynamic compliance.

Shaker

Ram tip acceleration

Figure 1

 1 Figure 1 General control structure of a rack and pinion drive

Figure 2

 2 Figure 2 Control scheme with Linear Fractional Transformation.



Figure 3 a

 3 Figure 3 a) vertical turning center, b) measured FRF, c) predicted FRF.

Figure 4 a

 4 Figure 4 a) FRF with active damping, b) Vibration amplitude during a cutting test.

Acknowledgment

The authors would like to acknowledge the financial support of the NSERC CANRIMT2 Canadian Project.


Q A Wang

L Nivanen

A El Kaabouchi

J P Badiali

A Le Méhauté

Incomplete normalization of probability on multifractals

Keywords: numbers : 05.20.-y (Classical statistical mechanics) 05.45.Df (Fractals) 02.50.-r (Probability theory, stochastic processes, and statistics)

This work is an extension of the incomplete probability theory from the simple case of monofractals previously studied to the more general case of multifractals which can occur in the phase space without equiprobable partition.

1) Introduction

The incomplete probability distribution (IPD) has been proposed for the case where the probability cannot be normalized, i.e., 1

1 ≠ ∑ = w i i p
. For a class of IPD, we have suggested the following incomplete normalization (IN) [1]:

1 1 = ∑ = w i q i p ( 1
)
The reader can find discussion of the various reasons for the existence of this nonadditive or unnormalizable probability distribution in the references [1][3][4] [START_REF] Nivanen | [END_REF][6] [7].

One of the cases for which IN is suggested is the chaotic systems evolving in phase space having fractal attractors [START_REF] Robert | Chaos and Nonlinear Dynamics[END_REF][3] [4]. In our previous papers [3][4], IN was studied on monofractals having homogeneous distribution of segments of same length such as the standard Cantor set in Figure 1. This set is constructed by iteration. At the k th iteration, the initial line of length L 0 is transformed into a set with N k = 2 k segments of length

0 3 1 L k k ki ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = δ = δ
(i = 1, 2, …, N k are in the order say from left to the right). To understand the link between this structure evolution and the probability distribution of nonequilibrium systems, we can consider an ensemble of chaotic systems which evolve in their phase space [START_REF] Beck | Thermodynamics of chaotic systems[END_REF] and are gradually attracted in a fractal structure (strange attractor) formed by their trajectories. At the same time, the occupied phase volume gradually maps into these attractors. To illustrate this evolution, let us take the example of a phase attractor in the form of the Cantor set of Figure 1. Suppose L 0 is the volume of the initial states, the total occupied volume at a later time t k of the k th iteration is

0 3 2 L N L k k k k ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = δ =
. If all the state points are equally probable, the usual definition of probability that a segment of length ki δ is visited by the system is given by

k ki ki L P δ =
. This probability can be normalized as follows :

1 1 1 = δ = δ = ∑ ∑ = = k k k k N i ki N i ki L N L P k k (2) Figure 1, Standard Cantor set (m = 2, ρ = 1/3) built at iteration orders from k = 0 to k = 4. At first iteration, the total length of the curve is 0 1 3 2 L L ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = with N 1 = 2 segments of length 0 1 3 1 L ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = δ ; at second iteration, 0 2 2 3 2 L L ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = with N 2 = 2 2 segments of length 0 2 2 3 1 L ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = δ , etc.
The probability defined in Eq.( 2) is scientifically sound and in agreement with the probability theory. Its only drawback is that it is not convenient for describing physical systems out of equilibrium. If the system is in equilibrium at a stage k, the above definition of probability is reasonable and sufficient for the statistical description of the equilibrium state on the total phase volume L k . However, if the system is out of equilibrium and L k is only an intermediate phase volume changing in time, it will be more convenient, in order to take into account the time evolution of the distribution, to define the probability finding a system in ki δ with respect to the

initial volume L 0 [1][3] , i.e., 0 L p ki ki δ = . ( 3
)
This probability not only gives us the nonequilibrium distribution at any moment but also the evolution of the distribution with respect to the initial conditions. Notice that this probability is not normalizable since we have 1

0 0 1 1 ≠ = δ = ∑ ∑ = = L L L p k N i ki N i ki k k .
In our previous work [START_REF] Wang | [END_REF], it was shown that, for the case of self-similar monofractal, the probability of Eq.( 3) can be normalized by :

1 1 = ∑ = k N i q ki p (4)
for any stage k with a unique q = d, where d is the self-similarity dimension of the monofractals (see definition below). The parameter q characterizes the evolution of the phase space. When q < 1, the phase volume shrinks in time, when q > 1, there is phase space expansion. The normalization becomes complete for q = 1 without evolution of distribution. Eq.( 4) has been derived on monofractals having segments of same length such as in Figure 1. It was only conjectured [START_REF] Wang | [END_REF] that the same calculation would hold for multifractals which in general contain segments of different lengths, but no proof has been given.

In this work, it will be shown that Eq.( 4) holds for any multifractal with self-similar structure having the same iteration rule at each stage. In order to be clear for the readers who have not followed the previous work, we will begin by the simplest case of homogeneous monofractals.

More complicated cases with different segments will be analyzed first with two segments of different lengths replacing a segment of precedent stage, and then with an arbitrary number m of different segments. Proof for the uniqueness of q will also be given.

2) Monofractals

In this section, we recapitulate some essential results obtained previously for the case of monofractals. Let us consider a deterministic mass fractal form F such as the usual Cantor set on Figure 1. Let F k be the intermediate k th stage of construction. When the initiator F 0 is a segment of length L 0 , F 1 is defined by N 1 = m copies of F 0 with same scale factor ρ. We suppose that two different copies never intersect. At the

k th stage, F k is composed of N k = m k segments of same length 0 L k k ki ρ = δ = δ
. The total length of the k th stage F k form is given by :

( ) 0 1 L m N L k k k N i ki k k ρ = δ = δ = ∑ = (5)
The notion of self-similar dimension is defined from the idea of fractal content of order α of the form F k :

( ) α α α α = α α ρ = ρ = δ = δ = ∑ 0 0 1 , ) ( L m L m N C k k k k k N i ki k k (6)
The fractal similarity dimension d is derived from the behavior of the content in the k → +∞ limit :

⎪ ⎩ ⎪ ⎨ ⎧ +∞ → < → > α α α α , , : 0 : k k C d C d (7)
In the monofractal case,

1 = ρ d m so that ) / 1 ln( ln ρ m d = .
Now let us consider two approaches leading to two different measures.

First we only look at the stage k and choose a point

M on F k . As F k is composed of N k segments, M
belongs to one of these segments. As mentioned above and used in many textbooks [START_REF] Beck | Thermodynamics of chaotic systems[END_REF], the measure

k ki ki L P δ =
is a probability. The monofractal character of F implies that at any stage

k k k ki P L P = δ = is independent of i.
It is given by :

k k k k ki k ki m L m L L P P 1 ) ( 0 0 = ρ ρ = δ = = (8)
and normalized by :

1 1 1 = = = ∑ = k k k k N i ki m m P N P k . ( 9
)
This distribution {P ki }, 1 ≤ i ≤ N k , only concerns the k th stage F k of the monofractal form. The fractal character of the distribution and of its time evolution is not taken into account. Now let us look at the fractal building in an evolutionary point of view by choosing a point M on the initiator F 0 and following its behavior in time. In a mass fractal like Cantor set, we have

F k ⊂ F k-1 for any k, which implies F k ⊂ F 0 .
So M may or may not belong to the k th stage F k . This is a kind of loss of state points (or information). To evaluate this evolution, we introduce the quantity 0 L p ki ki δ = [START_REF] Wang | [END_REF]. The monofractal character of F implies that at any stage

k k ki p L p = δ = 0 is independent of i.
We get from the definition of fractal dimension :

k kd d k d k d ki m L L L p 1 0 0 0 = ρ = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ρ = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ δ = (10)
and

1 1 1 = = = ∑ = k k d k k N i d ki m m p N p k (11)
The distribution {p ki }, 1 ≤ i ≤ N k is essentially different from the normalizable probability distribution Eq.( 8). As indicated above, it measures the probability to encounter the i th segment of F k when the system comes from somewhere on the initial set F 0 . It takes into account the fractal nature of the evolution of set F.

3) Binary multifractals

As in the case of monofractal sets, a probability law with IN can be defined on a multifractal set characterized by different scale factors for each copy of the initial segment. In this section, we first treat the case of a "binary" multifractal set. At each stage, a segment is replaced by m = 2 copies with scale factors ρ 1 and ρ 2 , respectively (0 < ρ 1 , ρ 2 < 1). An example is given on Figure 2.

(m = 2, ρ 1 = 1/2, ρ 2 = 1/3) built at iteration orders from k = 0 to k = 4. At the k th stage, F k is composed of N k = 2 k segments. Let j ∈ {0,…, k}. On F k , we can find ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ j k segments of length 0 2 1 , L j j k j k ρ ρ = Δ -
. The total length of the k th stage F k form is given by :

( ) 0 2 1 0 , 1 L j k L k k j j k N i ki k k ρ + ρ = Δ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = δ = ∑ ∑ = = (12)
The notion of self-similarity dimension can be extended to the case of a multifractal set. The content of order α of the form F k takes the form :

∑ ∑ ∑ = α α α - = α = α α ρ ρ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = Δ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = δ = k j j j k k j j k N i ki k L j k j k C k 0 0 2 ) ( 1 0 , 1 , (13) ( ) ( ) ( ) α α α = α α - α α ρ + ρ = ρ ρ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = ∑ 0 2 1 0 0 2 1 , L L j k C k k j j j k k (14)
The multifractal similarity dimension d is introduced by Eq.( 7). It verifies :

1 2 1 = ρ + ρ d d (15)
The quantities P ki and p ki are introduced in the same way as in Section 2, by setting respectively. We note that the distributions {P ki } and {p ki }, 1 ≤ i ≤ N k are no longer homogeneous due to the different scale factors ρ 1 and ρ 2 . Indeed, they verify the same normalization relationships as Eqs.( 9) and (11) in Section 2. Again, P ki is normalizable :

1 1 1 1 1 = = δ = ∑ ∑ = = k k N i ki k N i ki L L L P k k (16)
For the calculation of the sum of p ki d , let us classify all segments (denoted by i) of F k into groups (denoted by j) of segments of same length j k , Δ , such that the sum over i be replaced by a sum over

j, i.e., ∑ ∑ ∑ ∑ = - = = = ρ ρ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = Δ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = δ = k j d jd d j k d k j d j k d N i d ki d N i d ki L j k L j k L L p k k 0 0 2 ) ( 1 0 0 , 0 1 0 1 1 1 1 (17) or ( ) ( ) ( ) 1 1 2 1 0 2 1 1 = = ρ + ρ = ρ ρ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = ∑ ∑ = - = k k d d k j j d j k d N i d ki j k p k . ( 18
)
So in the binary multifractal case, a probability distribution with incomplete normalization such as Eq.( 1) can still be defined according to Eq.( 18) at each stage of a given multifractal with q=d.

Since the factor ( )

k 2 1 ρ + ρ
appears in Eq.( 12), all the distributions {P ki }, 1 ≤ i ≤ N k can be derived from the distribution of F 1 stage {P 1i }, 1 ≤ i ≤ 2 = N 1 . For instance, from Figure 2, we get respectively. In the same way, the factor ( )

P 11 = 3/
k d d 2 1 ρ + ρ
appears in Eq.( 18), so all the distributions

{p ki }, 1 ≤ i ≤ N k can be derived from the F 1 stage distribution {p 1i }, 1 ≤ i ≤ 2 = N 1 . From

4) General multifractals

In this section, we consider the case of general deterministic multifractal sets. At each stage, each segment is replaced by m copies with scale factors ρ 1 , ρ 2 , …, ρ m respectively. For any integer i ∈ {1, 2, …, m}, 0 < ρ i < 1. Such an example is given on Figure 3. At the k th stage, F k is composed of N k = m k segments of length δ ki . We note J any set of integers {j 1 , j 2 ,…, j m }, that verify :

{ } { } k j k j m i m i i i = ∈ ∈ ∀ ∑ =1 , ,..., 1 , 0 , ,..., 2 , 1 (19)
For any set J that verifies Eq.( 19), we can find

∏ = m i i j k 1 ! ! segments of length 0 1 , L m i j i J k i ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ρ = Δ ∏ = on F k .
The total length of the F k form at k th stage is given by :

0 1 1 1 ! ! L j k L J m i j i m i i N i ki k i k ∑ ∏ ∏ ∑ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ρ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ = δ = = = = (20)
where the sum is extended upon the all possible sets of the form in Eq.( 19). We see a multinomial series expansion such as

0 1 L L k m i i k ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ρ = ∑ = (21)
The content of order α of the form F k can be performed :

( ) α α α α α α α α ρ ρ ρ δ 0 1 0 1 1 0 1 1 1 , ! ! ! ! L L j k L j k C k m i i J m i j i m i i J m i j i m i i N i ki k i i k ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∑ = ∑ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∏ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ ∏ = ∑ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∏ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ ∏ = ∑ = = = = = = = (22)
We deduce the identity verified by the multifractal similarity dimension d :

1 1 = ρ ∑ = m i d i (23)
In order to see the uniqueness of d, let us define f by

∑ = α ρ = α m i i f 1 ) (
which is a strictly decreasing and continuous function. We have in addition f(α) = m for α=0 and f(α) → 0 for α → +∞, meaning that Eq.( 23) is satisfied by a unique real d.

The quantities P ki and p ki are defined in the same way as in Sections 2) and 3), by setting

k ki ki L P δ = and 0 L p ki ki δ =
respectively. The incomplete normalization relationship Eq.( 18) can be generalized into :

( ) ∑ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∏ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ ∏ = ∑ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∏ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ ∏ = ∑ Δ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ ∏ = ∑ = ∑ = = = = = = = J m i j d i m i i d J d m i j i m i i d J d J k m i i d d N i d ki N i d ki i i k k j k L j k L j k L L p 1 1 0 1 1 0 , 1 0 0 1 1 ! ! ! ! 1 ! ! 1 ρ ρ δ (24)
which is a multinomial series expansion :

1 1 1 1 = = ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ρ = ∑ ∑ = = k k m i d i N i d ki k p (25)
As in the binary case, an IPD with a unique q can be defined on a general deterministic multifractal set with the same rule of iteration at each stage. For all the cases considered up to now, the uniqueness of q comes from the uniqueness of the fractal dimension d. The multifractal case is particularly interesting because it makes it possible to extend incomplete distribution beyond the simple case of equiprobable phase space partition.

5) Multifractals with multidimensional initiators

The above cases concern only one dimensional initiator F 0 . In this section, we address the multifractals having multidimensional initiators. The multidimensional monofractals have been studied in a previous work [START_REF] Wang | [END_REF], which we will recapitulate in what follows before applying it to the multifractal case. An example of multidimensional monofractal defined by a unique scale factor ρ from a square initiator is given in Figure 4. The fractal dimension is still introduced from a length measure. In Figure 4, the length

0 L L L k k ki ρ = =
is given by the edge of any square of the fractal form F k . For a more general initiator, we choose a characteristic length and call it a gauge. The fractal dimension is still defined by Eqs.( 6) and ( 7).

In the plane, F 0 is a polygon of surface S 0 and

F k is composed of N k = m k polygons of surface 0 2 S k k ki ρ = σ = σ
. The total surface S k of the k th stage is given by :

( ) 0 2 1 S m N S k k k N i ki k k ρ = σ = σ = ∑ = (26)
The quantities P ki and p ki are now defined with respect to the surfaces, i.e., Let S 0 be the surface of the initiator F 0 . At the k th iteration, the form F k is composed of N k = m k copies of F 0 of surface ki σ . For any set J that verifies Eq.( 19), we can find

∏ = m i i j k 1 ! ! copies of F 0 of surface 0 1 2 , S S m i j i J k i ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ρ = ∏ = on F k .
The total surface of the k th stage F k form is given by : 0

1 2 1 1 ! ! S j k S J m i j i m i i N i ki k i k ∑ ∏ ∏ ∑ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ρ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ = σ = = = = (28)
where the sum is carried out over all the possible sets of the form Eq.( 19). We get a multinomial series expansion :

0 1 2 S S k m i i k ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ρ = ∑ = . ( 29
)
The incomplete normalization relationship Eq.(25) now becomes :

∑ ∏ ∑ ∑ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ = σ = = = = J d J k m i i d d N i d ki N i d ki S j k S S p k k 2 / , 1 2 / 0 2 / 0 1 2 / 1 2 / ! ! 1 (30) i.e., ( ) ∑ ∏ ∏ ∑ ∏ ∏ ∑ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ρ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ρ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ = = = = = = J m i j d i m i i d J d m i j i m i i d N i d ki i i k j k S j k S p 1 1 2 / 0 2 / 1 2 1 2 / 0 1 2 / ! ! ! ! 1 (31)
which is a multinomial series expansion : An incomplete probability distribution with unique q=d/2 (still resulting from the uniqueness of the fractal dimension) can be defined on a general bidimensional multifractal set with the same rule of iteration at each stage. This result can be extended to higher dimensional initiators. For the "three dimensional" case for example, one obtains an exponent q=d/3. For arbitrary topological dimension d T , we should have :

6) Concluding remarks

The present work is an extension of the IPD from the simple case of monofractals to the multifractals. This extension is useful for the understanding of the incomplete statistics and of their relevance in the study of nonequilibrium complex systems evolving in fractal geometry such as the strange attractors in phase space. The conclusion of the present work is that the mathematical formalism of IPD having incomplete normalization with a unique parameter q=d/d T is valid for any deterministic multifractal with the same rule at each iteration.

However, the uniqueness of q may be perturbed in the case of the fractals which change rules during their iterations or for different scales. The possible relationships between q and the different iteration rules and the consequences of this multi-rule evolution are under investigation. The result is expected to be helpful for the understanding of the evolution of many complex systems with obvious dynamical transition in time and in scale, such as the degree distributions in multiple regimes observed in many scale free networks.

Figure 2 ,

 2 Figure 2, Modified binary multifractal Cantor set (m = 2, ρ 1 = 1/2, ρ 2 = 1/3) built at iteration

  5, P 12 = 2/5 and P 21 = 9/25 = P 11 2 , P 22 = P 23 = 6/25 = P 11 P 12 , P 24 = 4/25 = P 12 2 ,

Figure 2 ,

 2 we have p 11 = 3/5, p 12 = 2/5, and p 21 = 9/25 = p 11 2 , p 22 = p 23 = 6/25 = p 11 p 12 , p 24 = 4/25 = p 12 2 . It means that for a given fractal F, the essential relationships for Eqs.(16) and (18) to hold are the iteration rules applied to k =1 and to each following stage of the structure.

Figure 3 ,

 3 Figure 3, Modified multifractal Cantor set (m = 3, ρ 1 = 1/3, ρ 2 = 1/6, ρ 3 = 1/4) built at iteration orders from k = 0 to k = 3.

Figure 4 ,

 4 Figure 4, bidimensional Cantor set (m = 2, ρ = 1/3) built at iteration orders from k = 0 to k = 3.

  distribution {p ki } still verifies an incomplete normalization relationship with a new exponent d/2.

Figure 5 ,

 5 Figure 5, bidimensional multifractal Cantor set (m = 4, ρ 1 = 1/5, ρ 2 = 1/3, ρ 3 = 1/2, ρ 4 = 1/4) built by iteration orders from k = 0 to k = 3.