# Integrated Multiobjective  Optimal Design for Active Control System Based on Genetic Algorithm

Authors: Ma Yong-Quan; Qiu Hong-Xing
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

The integrated multiobjective optimal design method for structural active control system is put forward based on improved Pareto multiobjective genetic algorithm, through which the position of actuator is synchronously optimized with active controller. External excitation is simulated by stationary filtered white noise. The root-mean-square (RMS) of structural response and active control force can be achieved by solving Lyapunov equation in the state space. The design of active controller adopts linear quadratic regulator (LQR) control algorithm. Minimum ratio of the maximum RMS of controlled structural displacement divided by the maximum RMS of uncontrolled structural displacement and minimum ratio of the maximum RMS of controlled structural shear divided by the maximum RMS of uncontrolled structural shear, together with minimization of the sum of RMS of active control force, are used as the three objective functions of multiobjective optimization. The optimization process takes the impact of structure and excitation parameter on the optimized results. An eight-storey six-span plane steel frame was used as an emulational example to demonstrate the validity of this optimization method. Results show that the proposed integrated multiobjective optimal design method is simple, efficient, and practical with good universality.


---

## Body

## 1. Introduction

In the research area of structural active control, optimal design of the control system, especially the position optimization of actuator, is always one of the hottest hot points of research. In nature, the position optimization of actuator is an issue of the optimization of discrete domain [1]. As a result, it is difficult to find the optimal solution through traditional methods [2]. Random search methods, such as simulated annealing algorithm [3, 4] and genetic algorithm [5, 6], are global optimization algorithms which can be used for the issues of discrete optimization. Thus, they are widely applied in the research of position optimization of control devices. However, previous researches have neglected the internal relations between the position optimization of control devices and the controller itself [7], and the optimization processes are always carried out in time domain [8, 9], in which the optimal position of actuator is worked out only based on certain external excitation and optimization rules [10, 11]. Due to the low computing efficiency, researchers tend to apply a simplified shearing series multiparticle model [12, 13]. On such condition, only the optimal storey of actuator instead of its detail plane position can be figured out [14], and the guiding significance of optimized results is inadequate in actual engineering application.Considering the high cost of initial investment and operation of structural active control system [15–17], the simplistic design may cause horrible waste of resources. In the opinion of the authors of this paper, the security, in the optimal design of active control system, improvement of functions, and economical efficiency of control system should be synthetically considered during its complete lifecycle. For this reason, the optimal design of active control system is essentially a matter of multiobjective optimization with universality. Combining improved Pareto optimal solution theory and genetic algorithm is an effective way for multiobjective optimization.This paper, based on improved Pareto multiobjective genetic algorithm, puts forward an active control integrated optimal design method. Compared with the method proposed in previous literatures, the method of optimization design presented in this paper, which can consider the internal relations between the control device’s position optimization and the controller itself, can obtain the optimal position of actuator under the condition of the random excitation and stochastic optimization criterion. The optimization process of the method proposed in this paper can be carried out in the frequency domain, so that the computational efficiency is very high, which is applicable to the complex plane model’s optimization of active control system, and can obtain the optimal layout position of actuator. In the meanwhile, this has important reference value for the practical application in engineering.In order to improve computing efficiency, the optimization process is carried out in frequency domain, the random seismic excitation is simulated by stationary filtered white noise. The variance of structural response and active control force can be worked out by obtaining the solution of Lyapunov equation of the structure under stationary filtered white noise in the state space. The design of active controller adopts linear quadratic regulator (LQR) controlling algorithm. Minimization of the ratio between the maximum root-mean-square (RMS) of controlled structural displacement and the maximum RMS of uncontrolled structural displacement and minimization of the ration between the maximum RMS of controlled structural shear and the maximum RMS of uncontrolled structural shear, together with minimization of the sum of RMS of active control force, are used as the three objective functions of multiobjective optimization. At the end, this paper takes a plane frame structure as an example to certify the effectiveness of proposed active control integrated multiobjective optimization design method.


## 2. Solving the Random Response Based on Lyapunov Equation

### 2.1. Variance of Uncontrolled Structural Response

The structural dynamic equation ofn-particle model excited by earthquake acceleration x¨g(t) is described as follows:
(1)Mx¨(t)+Cx˙(t)+Kx(t)=-Mμx¨g(t),
where M, C, and K are the structural mass matrix, damping matrix, and stiffness matrix, respectively; x(t), x˙(t), and x¨(t) are the structural displacement vector, velocity vector, and acceleration vector, respectively; μ is the unit column vector.The state equation of formula (1) can be written as
(2)Y˙=HY+Bx¨g,(3)Z=AY+Lx¨g,
where H=[0n×nIn×n-inv(M)×K-inv(M)×C], B=[0n×1-In×1], and Y=[xx˙] are structural state vectors; Z is the output vector; I and 0 are identity matrix and zero matrix, respectively, and the subscript is the dimensions corresponding to the matrix.The Clough-Penzien spectrum [18] of filtered white noise is used to simulate the input earthquake acceleration x¨g(t), and it can be presented by state equation as follows:
(4)Y˙s=HsYs+DsS0(5)Zs=AsS0,
where Ys=[xgY˙g], Hs=[01-ωg2-2ξgωg], Ds=[0-1], and As=[-ωg2-2ξgωg]T; ωg and ξg are the predominant frequency and damping ratio of foundation soil, respectively; S0 refers to the input white noise.The state equation of earthquake input-structural model can be obtained as follows by combining formulas (2) and (4) as well as formulas (3) and (5):
(6)Y˙τ=HτYτ+BτS0,(7)Zτ=AτYτ+LτS0,
where Hτ=[HBAs02×2nHs], Bτ=[02n×1Bs], As=[ALAs01×2nAs], Lτ=0(nz+1)×1, and the extended state vector Yτ=[YYs].The solution of differential equation (6) is
(8)Yτ(t)=exp⁡⁡(Hτ(t-t0))Yτ0+∫t0texp⁡(Hτ(t-T))BτS0(T)dT,
where Yτ0 is the initial condition when t=t0; exp⁡(Hτt) is the general solution of first-order equation; the integral term in the right side of the formula above is a particular solution of the equation. For convenience, given Yτ0=0, the formula above can be rewritten as
(9)Yτ(t)=∫t0texp⁡⁡(Hτ(t-T))BτS0(T)dT.The related functional matrix of structural response can be calculated through the formula above, and its expression is given as follows:(10)RYτYτ(t1,t2)=E[Yτ(t1)YτT(t2)]=∫t0t1∫t0t2exp⁡⁡(Hτ(t1-T1))BτE[S0(T1)S0τ(T2)]×BτTexp⁡(HτT(t2-T2))dT1dT2.In the formula above, the variance matrix of input excitation is defined as follows:(11)E[S0(T1)S0τ(T2)]=2πP0δ(t2-t1)=Psδ(t2-t1).If we plug the formula above into formula (10) and take the integral property of function δ into account, formula (10) can be transformed into
(12)RYτYτ(t1,t2)=∫t0t1exp⁡⁡(Hτ(t1-T1))BτPsBτTexp⁡(HτT(t2-T2))dT1.Provided thatt2=t1=t, the variance matrix BYτ of response is
(13)BYτ=∫t0texp⁡⁡(Hτ(t-T))BτPsBτTexp⁡⁡(HτT(t-T))dT=exp⁡⁡(Hτt)·(∫t0texp⁡⁡(-HτT)BτPsBτTexp⁡⁡(-HτTT)dT)exp⁡(HτTt).Taking the derivative of timet on both sides of the formula above, the result is
(14)dBYτdt=HτBYτ(t)+BYτ(t)HτT+BτPsBτT,
where formula (14) is called Lyapunov equation [19] that describes the property of response variance BYτ(t). The response maintains stability, and variance is irrelevant to time in a steady situation. At this point, formula (14) can be transformed into
(15)HτBYτ+BYτHτT=-Bτ2πP0(S)BτT.VarianceBYτ(t) of state vector Yτ can be obtained by solving the formula above and the variance of output vector Zτ is
(16)BZτ=AτBYτAτT+LτP0(S)LτT.


### 2.2. Variance of Both Controlled Structural Response and Control Force

The structural dynamic equation ofn-particle model installed with actuators excited by earthquake acceleration x¨g(t) can be expressed as
(17)Mx¨(t)+Cx˙(t)+Kx(t)=-Mμx¨g(t)+NUb.In this formula,N is the position matrix of control force; Ub refers to active control force. The following extended state equation can also be achieved according to the procedures in Section 2.1:
(18)Y˙τ=HτYτ+BτS0+DmUb,Zτ=AτYτ+LτS0+ΛτUb,Zm=AmYτ+LmS0+ΛmUb.Optimal feedback gain matrixQb is designed using LQR controller [20] so as to obtain expected optimal active control force Ub. And its quadratic objective function is expressed with the following form:
(19)V=∫0∞{ZYTGZY+UbTJUb}dt,
where ZY means control objective vector of the system; G and J are the weight matrixes weighing the structural response and active control force, respectively.The controller design and state estimation are treated separately in accordance with separation principle [21], and optimal control law is determined based on linear quadratic optimal control theory [22] as follows:
(20)Ub=-ΩbYb.In this formula,Ωb is the full state feedback gain matrix; Yb refers to the system state vector estimated based on Kalman filter [23], and it can be expressed as
(21)Y˙b=HτYb+DτUb+μb(Zm-AmYb-ΛmUb),
where μb is observer gain matrix and Zm means measured response. System state is estimated through measuring acceleration in this paper. Finally, the LQR controller can be described with the following state equation:
(22)Y˙b=HbYb+DbZm,Ub=AbYb-ΛbZm,
where Hb=Hτ-DτΩb-μbAm+μbΛmΩb, Db=μb, Ab=-Ωb, and Λb=0.Figure1 shows the block diagram of filtered white noise-controlled structure-LQR controller.Figure 1
Block diagram of stationary filtered white noise-controlled structure-controller.An extended state equation can be obtained through expanding controlled structure, LQR controller, and filtered white noise into a state equation, which is expressed as(23)Y˙d=HdYd+Dd[S0Ub],Zτ=AdYd-Λd[S0Ub].In this formula,Yd=[YTYsTYbT] refers to the system state; Hd, Db, Ad, and Λd represent the coefficient matrixes, and they can be obtained by using “feedback” command in MATLAB.Similarly, the variances of both the controlled structure response and the control force under stationary random excitation can be worked out. The peak control force needed in actual engineering application can be figured out conveniently through RMS and peak factor.


## 2.1. Variance of Uncontrolled Structural Response

The structural dynamic equation ofn-particle model excited by earthquake acceleration x¨g(t) is described as follows:
(1)Mx¨(t)+Cx˙(t)+Kx(t)=-Mμx¨g(t),
where M, C, and K are the structural mass matrix, damping matrix, and stiffness matrix, respectively; x(t), x˙(t), and x¨(t) are the structural displacement vector, velocity vector, and acceleration vector, respectively; μ is the unit column vector.The state equation of formula (1) can be written as
(2)Y˙=HY+Bx¨g,(3)Z=AY+Lx¨g,
where H=[0n×nIn×n-inv(M)×K-inv(M)×C], B=[0n×1-In×1], and Y=[xx˙] are structural state vectors; Z is the output vector; I and 0 are identity matrix and zero matrix, respectively, and the subscript is the dimensions corresponding to the matrix.The Clough-Penzien spectrum [18] of filtered white noise is used to simulate the input earthquake acceleration x¨g(t), and it can be presented by state equation as follows:
(4)Y˙s=HsYs+DsS0(5)Zs=AsS0,
where Ys=[xgY˙g], Hs=[01-ωg2-2ξgωg], Ds=[0-1], and As=[-ωg2-2ξgωg]T; ωg and ξg are the predominant frequency and damping ratio of foundation soil, respectively; S0 refers to the input white noise.The state equation of earthquake input-structural model can be obtained as follows by combining formulas (2) and (4) as well as formulas (3) and (5):
(6)Y˙τ=HτYτ+BτS0,(7)Zτ=AτYτ+LτS0,
where Hτ=[HBAs02×2nHs], Bτ=[02n×1Bs], As=[ALAs01×2nAs], Lτ=0(nz+1)×1, and the extended state vector Yτ=[YYs].The solution of differential equation (6) is
(8)Yτ(t)=exp⁡⁡(Hτ(t-t0))Yτ0+∫t0texp⁡(Hτ(t-T))BτS0(T)dT,
where Yτ0 is the initial condition when t=t0; exp⁡(Hτt) is the general solution of first-order equation; the integral term in the right side of the formula above is a particular solution of the equation. For convenience, given Yτ0=0, the formula above can be rewritten as
(9)Yτ(t)=∫t0texp⁡⁡(Hτ(t-T))BτS0(T)dT.The related functional matrix of structural response can be calculated through the formula above, and its expression is given as follows:(10)RYτYτ(t1,t2)=E[Yτ(t1)YτT(t2)]=∫t0t1∫t0t2exp⁡⁡(Hτ(t1-T1))BτE[S0(T1)S0τ(T2)]×BτTexp⁡(HτT(t2-T2))dT1dT2.In the formula above, the variance matrix of input excitation is defined as follows:(11)E[S0(T1)S0τ(T2)]=2πP0δ(t2-t1)=Psδ(t2-t1).If we plug the formula above into formula (10) and take the integral property of function δ into account, formula (10) can be transformed into
(12)RYτYτ(t1,t2)=∫t0t1exp⁡⁡(Hτ(t1-T1))BτPsBτTexp⁡(HτT(t2-T2))dT1.Provided thatt2=t1=t, the variance matrix BYτ of response is
(13)BYτ=∫t0texp⁡⁡(Hτ(t-T))BτPsBτTexp⁡⁡(HτT(t-T))dT=exp⁡⁡(Hτt)·(∫t0texp⁡⁡(-HτT)BτPsBτTexp⁡⁡(-HτTT)dT)exp⁡(HτTt).Taking the derivative of timet on both sides of the formula above, the result is
(14)dBYτdt=HτBYτ(t)+BYτ(t)HτT+BτPsBτT,
where formula (14) is called Lyapunov equation [19] that describes the property of response variance BYτ(t). The response maintains stability, and variance is irrelevant to time in a steady situation. At this point, formula (14) can be transformed into
(15)HτBYτ+BYτHτT=-Bτ2πP0(S)BτT.VarianceBYτ(t) of state vector Yτ can be obtained by solving the formula above and the variance of output vector Zτ is
(16)BZτ=AτBYτAτT+LτP0(S)LτT.


## 2.2. Variance of Both Controlled Structural Response and Control Force

The structural dynamic equation ofn-particle model installed with actuators excited by earthquake acceleration x¨g(t) can be expressed as
(17)Mx¨(t)+Cx˙(t)+Kx(t)=-Mμx¨g(t)+NUb.In this formula,N is the position matrix of control force; Ub refers to active control force. The following extended state equation can also be achieved according to the procedures in Section 2.1:
(18)Y˙τ=HτYτ+BτS0+DmUb,Zτ=AτYτ+LτS0+ΛτUb,Zm=AmYτ+LmS0+ΛmUb.Optimal feedback gain matrixQb is designed using LQR controller [20] so as to obtain expected optimal active control force Ub. And its quadratic objective function is expressed with the following form:
(19)V=∫0∞{ZYTGZY+UbTJUb}dt,
where ZY means control objective vector of the system; G and J are the weight matrixes weighing the structural response and active control force, respectively.The controller design and state estimation are treated separately in accordance with separation principle [21], and optimal control law is determined based on linear quadratic optimal control theory [22] as follows:
(20)Ub=-ΩbYb.In this formula,Ωb is the full state feedback gain matrix; Yb refers to the system state vector estimated based on Kalman filter [23], and it can be expressed as
(21)Y˙b=HτYb+DτUb+μb(Zm-AmYb-ΛmUb),
where μb is observer gain matrix and Zm means measured response. System state is estimated through measuring acceleration in this paper. Finally, the LQR controller can be described with the following state equation:
(22)Y˙b=HbYb+DbZm,Ub=AbYb-ΛbZm,
where Hb=Hτ-DτΩb-μbAm+μbΛmΩb, Db=μb, Ab=-Ωb, and Λb=0.Figure1 shows the block diagram of filtered white noise-controlled structure-LQR controller.Figure 1
Block diagram of stationary filtered white noise-controlled structure-controller.An extended state equation can be obtained through expanding controlled structure, LQR controller, and filtered white noise into a state equation, which is expressed as(23)Y˙d=HdYd+Dd[S0Ub],Zτ=AdYd-Λd[S0Ub].In this formula,Yd=[YTYsTYbT] refers to the system state; Hd, Db, Ad, and Λd represent the coefficient matrixes, and they can be obtained by using “feedback” command in MATLAB.Similarly, the variances of both the controlled structure response and the control force under stationary random excitation can be worked out. The peak control force needed in actual engineering application can be figured out conveniently through RMS and peak factor.


## 3. The Definition of Optimization Principle

Considering that the maximum output force of actuator is limited by its physical conditions in actual engineering application [24], the paper makes an assumption that the RMS control force meets the following constraints:
(24)max⁡{WUactivei}≤WUpermit.In this formula,WUpermit is the allowed maximum RMS of active control force, and WUactivei is the RMS of active control force generated by the ith actuator.The previous constraints can be treated with penalty function method, and the definition of optimized objective function is(25)∇1=max⁡j⁡{Wxjk}max⁡j⁡{Wxjwk}+η×max⁡[0,max⁡{WUpermitj}WUpermit-1],∇2=max⁡j⁡{WQjk}max⁡j⁡{WQjwk}+η×max⁡[0,max⁡{WUpermitj}WUpermit-1],∇3=∑j=1mWUactivej+η×max⁡[0,max⁡{WUpermitj}WUpermit-1].
In the formula, Wxjk and WQjk are the RMS of both lateral relative displacement and base shear on the jth storey of controlled structure, respectively; Wxjwk and WQjwk are the RMS of both lateral relative displacement and base shear on the jth storey of uncontrolled structure, respectively;m refers to the number of actuators; ∇1 represents the ratio between maximum RMS of lateral relative displacements of both controlled structure and uncontrolled structure; ∇2 means the ratio between maximum RMS of base shears of both controlled structure and uncontrolled structure; ∇3 refers to the sum of RMS of active control force. Constraints are taken into account in the second item; η represents penalty factor which is always expressed as a large constant. In the event that the constraints are satisfied, the value of this item is zero. Contrarily, the value of objective function is big while the possibility that it is chosen into next generation is small so as to inflict punishment on this individual. During the optimization process, the closed-loop stability of the system should be guaranteed; otherwise ∇1, ∇2, and ∇3 will be assigned with bigger values as a punishment to the individual.


## 4. Improved Multiobjective Genetic Algorithm Based on Pareto Theory

In this paper, improved multiobjective genetic algorithm based on Pareto optimal solution [25] is adopted to make multiobjective integrated optimal design of active control system. In order to obtain a better understanding of this algorithm, three definitions are introduced as follows.


### 4.1. Pareto Optimal Solution

The issues of multiobjective optimization above can be expressed as minimized objective function vector, which is expressed as(26)min⁡∇(y)={∇1(y),∇2(y),∇3(y)}.In this formula,y is the decision variable (refer to the positions of actuator and the design parameters of controller), ℧ is the decision variable space, Ω refers to objective function space, and Ω∈R3.Define optimal solutiony*∈℧. And given that there is at least one r∈X where X={1,2,3},
(27)∇r(y)>∇r(y*),
where y*∈℧ is the Pareto optimal solution. Its manifestation in objective function space is Pareto optimal boundary (known as P optimal frontier curve). As illustrated in Figure 2, Pareto optimal boundary is the curve made up by individuals whose Rank=1.Figure 2
Pareto ordering of two objectives.


### 4.2. Dominance Relation

Given thatλ and ϕ are any two different individuals of the evolving group, the two following conditions should be satisfied if λ dominates ϕ:(i)
for all the subobjectives,λ is not inferior to ϕ; namely, ∇r(λ)≤∇r(ϕ) is always true when ∀r∈X;(ii)
there should be at least one subobjective that makesλ superior to ϕ. Namely, ∇h(λ)<∇h(ϕ) is true when ∃h∈{1,2,3}.According to the definition of domination, it is easy to find out that the individuals on optimal boundary are not dominated. Provided that its boundary set sequence number is 1 (Rank=1), the boundary set sequence numbers of other individuals are defined as the numbers of individuals dominating it plus 1. As illustrated in Figure 2, the two subobjectives ∇1 and ∇3 form Pareto optimal curve and the three objectives make up an optimal surface. This curve divides objective function space into two parts, the upper space of which is feasible region and the lower space is infeasible region. As illustrated in Figure 2, for the curve coordinate (x0, y0), when the control effect to lateral displacement is 1-x0, the corresponding minimum (optimal) control force is y0, and such control effect cannot be realized if the force is lower than minimum control force.


### 4.3. Crowding Distance

In the event that sequence numbers of boundary set are equal, the pros and cons of individuals are represented by clustering distance. As illustrated in Figure2, for two subobjectives ∇1 and ∇3, the clustering distance of individual i is the sum of length and width of the dashed quadrangle in the figure. The crowding distance L~(i) of individual i can be defined as
(28)L~(i)=(Q~(i+1)·∇1-Q~(i-1)·∇1)+(Q~(i+1)·∇2-Q~(i-1)·∇2),
where Q~(i)·∇r represents the function value of individual i on the rth subobjective. When there are k subobjectives, (28) can be reexpressed as
(29)L~(i)=∑r=1k(Q~(i+1)·∇r-Q~(i-1)·∇r).Reference [26] puts forward multiobjective optimal algorithm (NSGA-II); however, this optimal algorithm NSGA-II, when facing the complex optimization problems, has relatively poor computational efficiency and tends towards premature convergence. To overcome the above shortcomings of NSGA-II, an improved multiobjective optimal algorithm is presented in this paper, and the partial ordering relation is defined as follows: during the generation process of new groups, if there are two individuals belonging to different boundary set, the individual with smaller boundary set sequence number is set as priority; if the two boundary sequence numbers are equal, the individual with larger clustering distance is set as priority. The flow diagram of this improved algorithm is shown in Figure 3, and the detailed descriptions are introduced as follows.(1)
Randomly generate initial parent populationD0 and progeny population F0, the population sizes of D0 and F0 are both N-, and define initial time t=0.(2)
Integrate the two populations (Dt and Ft) into St; that is, St=Dt∪Ft, and in accordance with the concept of Pareto construct nondomination level surface H of the St; that is, H={H1,H2,…,Hi}.(3)
EmptyDt+1 (Dt+1 is the parent population at the moment t+1) and set i=1.(4)
If|Dt+1|+|Hi|≤N-, then proceed to Step (5); otherwise, skip to Step (7).(5)
Calculate the clustering distance of individuals inHi, and add individuals in Hi to the parent population Dt+1; that is, Dt+1=Dt+1∪Ht.(6)
Seti=i+1; return to Step (4).(7)
Perform roulette selection, intersection, and variation operation onDt+1 to generate new population Ft+1.(8)
Ift<imax⁡ (imax⁡ is the maximum evolutionary number of generations), then return to Step (2); otherwise, the algorithm terminates.Figure 3
Flow diagram of improved NSGA-II algorithm.


## 4.1. Pareto Optimal Solution

The issues of multiobjective optimization above can be expressed as minimized objective function vector, which is expressed as(26)min⁡∇(y)={∇1(y),∇2(y),∇3(y)}.In this formula,y is the decision variable (refer to the positions of actuator and the design parameters of controller), ℧ is the decision variable space, Ω refers to objective function space, and Ω∈R3.Define optimal solutiony*∈℧. And given that there is at least one r∈X where X={1,2,3},
(27)∇r(y)>∇r(y*),
where y*∈℧ is the Pareto optimal solution. Its manifestation in objective function space is Pareto optimal boundary (known as P optimal frontier curve). As illustrated in Figure 2, Pareto optimal boundary is the curve made up by individuals whose Rank=1.Figure 2
Pareto ordering of two objectives.


## 4.2. Dominance Relation

Given thatλ and ϕ are any two different individuals of the evolving group, the two following conditions should be satisfied if λ dominates ϕ:(i)
for all the subobjectives,λ is not inferior to ϕ; namely, ∇r(λ)≤∇r(ϕ) is always true when ∀r∈X;(ii)
there should be at least one subobjective that makesλ superior to ϕ. Namely, ∇h(λ)<∇h(ϕ) is true when ∃h∈{1,2,3}.According to the definition of domination, it is easy to find out that the individuals on optimal boundary are not dominated. Provided that its boundary set sequence number is 1 (Rank=1), the boundary set sequence numbers of other individuals are defined as the numbers of individuals dominating it plus 1. As illustrated in Figure 2, the two subobjectives ∇1 and ∇3 form Pareto optimal curve and the three objectives make up an optimal surface. This curve divides objective function space into two parts, the upper space of which is feasible region and the lower space is infeasible region. As illustrated in Figure 2, for the curve coordinate (x0, y0), when the control effect to lateral displacement is 1-x0, the corresponding minimum (optimal) control force is y0, and such control effect cannot be realized if the force is lower than minimum control force.


## 4.3. Crowding Distance

In the event that sequence numbers of boundary set are equal, the pros and cons of individuals are represented by clustering distance. As illustrated in Figure2, for two subobjectives ∇1 and ∇3, the clustering distance of individual i is the sum of length and width of the dashed quadrangle in the figure. The crowding distance L~(i) of individual i can be defined as
(28)L~(i)=(Q~(i+1)·∇1-Q~(i-1)·∇1)+(Q~(i+1)·∇2-Q~(i-1)·∇2),
where Q~(i)·∇r represents the function value of individual i on the rth subobjective. When there are k subobjectives, (28) can be reexpressed as
(29)L~(i)=∑r=1k(Q~(i+1)·∇r-Q~(i-1)·∇r).Reference [26] puts forward multiobjective optimal algorithm (NSGA-II); however, this optimal algorithm NSGA-II, when facing the complex optimization problems, has relatively poor computational efficiency and tends towards premature convergence. To overcome the above shortcomings of NSGA-II, an improved multiobjective optimal algorithm is presented in this paper, and the partial ordering relation is defined as follows: during the generation process of new groups, if there are two individuals belonging to different boundary set, the individual with smaller boundary set sequence number is set as priority; if the two boundary sequence numbers are equal, the individual with larger clustering distance is set as priority. The flow diagram of this improved algorithm is shown in Figure 3, and the detailed descriptions are introduced as follows.(1)
Randomly generate initial parent populationD0 and progeny population F0, the population sizes of D0 and F0 are both N-, and define initial time t=0.(2)
Integrate the two populations (Dt and Ft) into St; that is, St=Dt∪Ft, and in accordance with the concept of Pareto construct nondomination level surface H of the St; that is, H={H1,H2,…,Hi}.(3)
EmptyDt+1 (Dt+1 is the parent population at the moment t+1) and set i=1.(4)
If|Dt+1|+|Hi|≤N-, then proceed to Step (5); otherwise, skip to Step (7).(5)
Calculate the clustering distance of individuals inHi, and add individuals in Hi to the parent population Dt+1; that is, Dt+1=Dt+1∪Ht.(6)
Seti=i+1; return to Step (4).(7)
Perform roulette selection, intersection, and variation operation onDt+1 to generate new population Ft+1.(8)
Ift<imax⁡ (imax⁡ is the maximum evolutionary number of generations), then return to Step (2); otherwise, the algorithm terminates.Figure 3
Flow diagram of improved NSGA-II algorithm.


## 5. Simulation Analysis

This paper takes an eight-layer six-span frame structure as an example to examine the effectiveness of new integrated multiobjective optimal algorithm for active control system. As illustrated in Figure4, numbers from 1~48 are possible positions of the control devices (actuators). Every node of this structure has three degrees of freedom, namely, 2 translational motions and 1 rotation. Both the elasticity modulus of the beam and that of the column are E=210 GPa, the column density ρc=7900 kg·m−3, and the beam density ρb=21000 kg·m−3. The piezoelectric ceramic actuators (PCA) are used as the active control devices in this paper. The nominal maximum output force of PCA is 1200 kN, and the nominal maximum stroke of PCA is 60 cm. See Table 1 for design parameters of beam and column.Table 1
Structural parameters.
Storey number
Column
Beam
Bending inertia momentI/m4
Cross-sectional areaA/m2
Bending inertia momentI/m4
Cross-sectional areaA/m2
1-2
1.10
×
10
-
4
2.20
×
10
-
2
8.10
×
10
-
4
2.20
×
10
-
2
3-4
1.00
×
10
-
4
1.90
×
10
-
2
7.60
×
10
-
4
1.90
×
10
-
2
5-6
0.85
×
10
-
4
1.60
×
10
-
2
6.80
×
10
-
4
1.60
×
10
-
2
7-8
0.75
×
10
-
4
1.40
×
10
-
2
6.00
×
10
-
4
1.40
×
10
-
2Figure 4
Finite element model of the structure installed with actuators.A specific set of reference values of parametersωg and ξg is given in [27] according to the site classification set forth in China. See Table 2 for parameters of four sorts of sites in the first group of grouping of design earthquake selected in this paper. The measured vector Zm takes the horizontal acceleration of each node as its value, and controlled vector ZY takes the horizontal displacement and velocity of each node as its value. The allowed maximum RMS of control force WUpermit=1.2×105 kN, and the weight matrixes G=10ϕ×I60×60 and J=Im×m, where m represents the number of actuators whose value is 10. The initial population size is set as 600, and the evolutionary number of generations N-s=60. Real coding is applied whose length is 49 dimensions, the first 48 of which are used to identify the position of actuators and the last one is used to identify the weight coefficient ϕ. In simulation analysis, the stiffness ratios of beam to column and site soil feature parameters have different effects on the optimization results of control system.Table 2
Design value of site soil parameters.
Parameter
Site classification
I
II
III
IV
ω
g
25.12
17.96
13.95
9.68
ξ
g
0.65
0.71
0.82
0.91Optimal frontier curve is illustrated in Figure5. For intuition and clearance of marking, only 50 individuals are depicted in a two-dimensional view, in which we can observe that the distribution of Pareto solution set is relatively even. More importantly, improved NSGA-II algorithm is able to provide a Pareto solution set compared to the fact that only one optimal solution can be determined with traditional optimization methods. This algorithm allows decision-makers to make flexible selections of layout schemes according to actual engineering situations, which manifests the superiority of improved active control multiobjective optimal algorithm put forward in this paper.Optimal frontier curve of Pareto.
(a)
∇
1
~
∇
3
(b)
∇
2
~
∇
3The actuators’ positions and relevant control gains corresponding to some optimal individuals selected from Pareto optimal frontier curve are listed in Table3 to Table 5, and they correspond to different stiffness ratios of beam to column and the site classifications. By comparing the optimized results shown in Tables 3 and 4, we will discover that site classifications have little impact on the optimized results of actuators’ positions if there is little difference between seismic reduction effect of displacement and seismic reduction effect of shear (1-∇1 and 1-∇2, resp.), which means the genetic algorithm has strong robustness. As illustrated in Figure 5 and Table 3 to Table 5, when stiffness ratio of beam to column kbc=0.45, the sum of RMS of control force (∇3) is a little smaller than that when kbc=4.5, and the ∇3 in site classification I is a little larger than that in site classification III. This is mainly caused by the changes of both structural dynamic characteristics and external excitation parameters. The sum of control force increases gradually as the control effect improves. Serial number 1 in Table 3 to Table 5 also indicates that actuators should be deployed on every storey so as to obtain good control effect and that most devices need to be deployed on the two bottom storeys. When the preset control effect declines to a certain range (which corresponds to Serial numbers 2 and three in Table 3 to Table 5), the optimal position of the actuator remains almost still, and active control effect is connected only with control gain, which can be achieved by adjusting weight matrix. The fourth row (Serial number 4) in Table 3 to Table 5 shows the results when actuators are deployed on three bottom storeys according to experience on condition that the weight matrix is the same as that in Serial number 1, which indicates that the system performance index is improved after optimization.Table 3
Optimized results of site classification I when linear stiffness ratio of beam to column is 4.5.
Serial number
Position of the actuator
∇
1
∇
2
∇
3/kN
ϕ
1
1, 2, 6, 7, 9, 12, 27, 33, 39, 43, 45, 48
0.23
0.37
481
10
9.89
2
7, 12, 25, 36, 37, 42, 43, 44, 45, 46, 47, 48
0.48
0.54
193
10
8.48
3
4, 12, 26, 35, 38, 41, 43, 44, 45, 46, 47, 48
0.69
0.73
92
10
7.75
4
1, 2, 4, 5, 6, 7, 8, 9, 11, 12, 13, 15
0.26
0.36
513
10
9.89Table 4
Optimized results of site classification III when linear stiffness ratio of beam to column is 4.5.
Serial number
Position of the actuator
∇
1
∇
2
∇
3/kN
ϕ
1
3, 5, 7, 10, 21, 29, 33, 37, 41, 44, 46, 48
0.22
0.34
405
10
10
2
12, 24, 34, 36, 38, 40, 42, 44, 45, 46, 47, 48
0.45
0.47
189
10
8.61
3
11, 24, 35, 36, 37, 39, 42, 44, 45, 46, 47, 48
0.63
0.65
112
10
7.96
4
1, 2, 4, 5, 6, 7, 8, 9, 11, 12, 13, 15
0.22
0.35
465
10
10Table 5
Optimized results of site classification III when linear stiffness ratio of beam to column is 0.45.
Serial number
Position of the actuator
∇
1
∇
2
∇
3/kN
ϕ
1
1, 2, 6, 7, 12, 27, 31, 36, 39, 43, 45, 48
0.27
0.35
312
10
9.27
2
12, 24, 34, 36, 38, 40, 42, 44, 45, 46, 47, 48
0.45
0.52
155
10
8.19
3
11, 24, 33, 36, 38, 39, 42, 44, 45, 46, 47, 48
0.64
0.70
82
10
7.49
4
1, 2, 4, 5, 6, 7, 8, 9, 11, 12, 13, 15
0.28
0.37
339
10
9.27It turns out by the simulation analysis that the optimal controller of active control system is closely related to the position of actuator and that they both are influenced by structural parameters and excitation characteristics. That is why the design of traditional active control system could not obtain optimal seismic reduction effect, because the optimization is done separately into system design optimization and action controller position optimization. As a result, optimal control effect can be obtained only when an integrated design in structure-control system is carried out.


## 6. Conclusions

Integrated multiobjective optimal design for active control system is presented in this paper by using Pareto optimal solution-based improved multiobjective genetic algorithm and random vibration theory. The conclusion is given as follows.(1)
It is necessary for active control system to make integrated optimal design of the control device position and the controller. Previous optimal design is done separately into controller design optimization and action controller position optimization, which is the major reason why optimal seismic reduction effect cannot be achieved.(2)
The Pareto optimal solution set based on improved Pareto multiobjective genetic algorithm provides decision-makers with even wider choices, with which they may select corresponding layout schemes in accordance with actual projects. Therefore, it has high project application value.With simplicity, efficiency, and practicability, integrated multiobjective optimization method for active control system put forward in this paper is expected to have promising engineering application.



---

*Source: 748237-2014-07-14.xml*