Loading [MathJax]/jax/output/SVG/jax.js
Research article

Blow-up properties of solutions to a class of p-Kirchhoff evolution equations

  • This paper is devoted to an initial-boundary value problem for a class of p-Kirchhoff type parabolic equations. Firstly, we consider this problem with a general nonlocal coefficient M(upp) and a general nonlinearity k(t)f(u). A new finite time blow-up criterion is established, also, the upper and lower bounds for the blow-up time are derived. Secondly, we deal with the case that M(upp)=a+bupp, k(t)1 and f(u)=|u|q1u, which was considered by Li and Han [Math. Model. Anal. 2019; 24: 195-217] only for q>2p1. The threshold results for the existence of global and finite time blow-up solutions to this problem are obtained for the case 1<q2p1, which, together with the results given by Li and Han, shows that q=2p1 is critical for the existence of finite time blow-up solutions to this problem. These results partially generalize and extend some recent ones in previous literature.

    Citation: Hui Yang, Futao Ma, Wenjie Gao, Yuzhu Han. Blow-up properties of solutions to a class of p-Kirchhoff evolution equations[J]. Electronic Research Archive, 2022, 30(7): 2663-2680. doi: 10.3934/era.2022136

    Related Papers:

    [1] Toufik Bakir, Bernard Bonnard, Jérémy Rouot . A case study of optimal input-output system with sampled-data control: Ding et al. force and fatigue muscular control model. Networks and Heterogeneous Media, 2019, 14(1): 79-100. doi: 10.3934/nhm.2019005
    [2] Paola Goatin, Chiara Daini, Maria Laura Delle Monache, Antonella Ferrara . Interacting moving bottlenecks in traffic flow. Networks and Heterogeneous Media, 2023, 18(2): 930-945. doi: 10.3934/nhm.2023040
    [3] Michael Herty, Lorenzo Pareschi, Sonja Steffensen . Mean--field control and Riccati equations. Networks and Heterogeneous Media, 2015, 10(3): 699-715. doi: 10.3934/nhm.2015.10.699
    [4] Antoine Tordeux, Claudia Totzeck . Multi-scale description of pedestrian collective dynamics with port-Hamiltonian systems. Networks and Heterogeneous Media, 2023, 18(2): 906-929. doi: 10.3934/nhm.2023039
    [5] Oliver Kolb, Simone Göttlich, Paola Goatin . Capacity drop and traffic control for a second order traffic model. Networks and Heterogeneous Media, 2017, 12(4): 663-681. doi: 10.3934/nhm.2017027
    [6] Dengfeng Sun, Issam S. Strub, Alexandre M. Bayen . Comparison of the performance of four Eulerian network flow models for strategic air traffic management. Networks and Heterogeneous Media, 2007, 2(4): 569-595. doi: 10.3934/nhm.2007.2.569
    [7] Rudy R. Negenborn, Peter-Jules van Overloop, Tamás Keviczky, Bart De Schutter . Distributed model predictive control of irrigation canals. Networks and Heterogeneous Media, 2009, 4(2): 359-380. doi: 10.3934/nhm.2009.4.359
    [8] Zhong-Jie Han, Gen-Qi Xu . Dynamical behavior of networks of non-uniform Timoshenko beams system with boundary time-delay inputs. Networks and Heterogeneous Media, 2011, 6(2): 297-327. doi: 10.3934/nhm.2011.6.297
    [9] Mattia Bongini, Massimo Fornasier . Sparse stabilization of dynamical systems driven by attraction and avoidance forces. Networks and Heterogeneous Media, 2014, 9(1): 1-31. doi: 10.3934/nhm.2014.9.1
    [10] Vincent Renault, Michèle Thieullen, Emmanuel Trélat . Optimal control of infinite-dimensional piecewise deterministic Markov processes and application to the control of neuronal dynamics via Optogenetics. Networks and Heterogeneous Media, 2017, 12(3): 417-459. doi: 10.3934/nhm.2017019
  • This paper is devoted to an initial-boundary value problem for a class of p-Kirchhoff type parabolic equations. Firstly, we consider this problem with a general nonlocal coefficient M(upp) and a general nonlinearity k(t)f(u). A new finite time blow-up criterion is established, also, the upper and lower bounds for the blow-up time are derived. Secondly, we deal with the case that M(upp)=a+bupp, k(t)1 and f(u)=|u|q1u, which was considered by Li and Han [Math. Model. Anal. 2019; 24: 195-217] only for q>2p1. The threshold results for the existence of global and finite time blow-up solutions to this problem are obtained for the case 1<q2p1, which, together with the results given by Li and Han, shows that q=2p1 is critical for the existence of finite time blow-up solutions to this problem. These results partially generalize and extend some recent ones in previous literature.



    Functional Electrical Stimulation (FES) consists of applying an electrical stimulation to the muscle, in order to produce functional movements. It can be used for the muscular reinforcement, reeducation of the muscle and in the case of paralysis to activate the paralyzed muscles. Mathematically FES leads to a sampled-data control problem which can be analyzed in this framework.

    The simulations of muscular response to electrical stimulations are based on dynamics models. The origin comes from the Hill-Langmuir equation in the context of biochemistry and pharmacology, see [12]. More recent models in the framework of model identification in non-linear control are due to Bobet and Stein [2], Law and Shields [15] and a more sophisticated model was proposed by Ding et al. [6,7,8,9] where the force model is coupled to a fatigue model based on experimental results [20]. This led to a set of five differential equations with a sampled-data control which can be used to describe the force response to a pulses train of electrical stimulations and we shall refer to this model as the Ding et al. force-fatigue model in the sequel.

    There is a limited literature that used this model to design an optimized train of pulses to control the force level [6,7,8,9,1]. Our aim is to produce a more complete study of the problem in the framework of non-linear optimal control, using sampled-data controls.

    A geometric analysis of the model is provided and preliminary results are presented in the framework of optimal control with sampled-data control [3]. It is based on a simplified dynamics using the geometric properties of the force-fatigue model and control reduction to simplify the physical control constraints. The complete system is analyzed in details using a control predictive strategy (MPC), see [17,19] coupled with a non-linear observer based on [11] to estimate the state variables.

    The article is organized as follows. In Section 2, we make a brief presentation of Ding et al. force-fatigue model based on [20]. In Section 3, the dynamics of the force model is investigated to describe the input-output properties. In Section 4, the force-fatigue model is analyzed in the framework of geometric optimal sampled-data control systems and preliminary results are presented with a simplified model using a model reduction and an input transformation. Section 5 is devoted to the observer description. In Section 6, MPC method is presented using a further discretization of the dynamics to conclude with the algorithm implemented to compute in practice the optimized pulses trains. Numerical results are presented in the final Section 7.

    We refer to [20] for a complete description and discussions of the model. The Ding et al. model studied in this article is presented next in the framework of model dynamics construction based on the so-called tetania phenomenon in muscular responses. The first part of the model is the force response (output) to electrical stimulations pulses (input). The pulses are normalized Dirac impulses at times 0=t1<t2<<tn:

    v(t)=ni=1δ(tti)

    and Ii=titi1 is the interpulse and convexifying leads to apply the input:

    v(t)=ni=1ηiδ(tti).

    with parameters ηi[0,1],i{1,,n}. Such pulses train feeds a first-order model to produce an output according to the dynamics

    ˙u(t)+u(t)τc=1τcni=1Riηiδ(tti) (1)

    where Ri is a scaling factor associated to the phenomenon of tetania and corresponds to an accumulated effect of successive pulses and is modelled as

    Ri={1 for i=11+(R01)exp(titi1τc) for i>1

    where the magnitude is characterized by R0 which is the limit case to high-frequency pulse. The FES response to the input is denoted by Es(t) and is given by

    Es(t)=1τcni=1RiηiH(tti)exp(ttiτc) (2)

    where H(tti)={0 if t<ti1 if tti is the Heaviside function.

    The force response to such a train is modelled by the two equations of the so-called force model:

    dCNdt(t)+CN(t)τc=Es(t) (3)

    which corresponds to a first order (resonant) linear dynamics which can be integrated with CN(0)=0 and the force response F(t) is described by the equation

    dFdt(t)=A(t)a(t)F(t)b(t) (4)

    where A is a fatigue variable. Non-linear features of the model are described by the two mappings a,b:

    a(t)=CN(t)Km(t)+CN(t),b(t)=1τ1(t)+τ2a(t) (5)

    where Km,τ1 are fatigue variables and τ2 an additional parameter.

    The complete model is obtained by describing the evolutions of the variables associated to fatigue and corresponds to the linear dynamics

    dAdt(t)=A(t)Arestτfat+αAF(t) (6)
    dKmdt(t)=Km(t)Km,restτfat+αKmF(t) (7)
    dτ1dt(t)=τ1(t)τ1,restτfat+ατ1F(t). (8)

    The full set of equations (3)-(4)-(6)-(7)-(8) are the force and fatigue model.

    We refer to Table 1 for the definitions and details of the symbols of the force-fatigue model.

    Table 1. 

    Margin settings

    .
    Symbol Unit Value description
    CN Normalized amount of
    Ca2+-troponin complex
    F N Force generated by muscle
    ti ms Time of the ith pulse
    n Total number of
    the pulses before time t
    i Stimulation pulse index
    τc ms 20 Time constant that commands
    the rise and the decay of CN
    R0 1.143 Term of the enhancement
    in CN from successive stimuli
    A Nms Scaling factor for the force and
    the shortening velocity
    of muscle
    τ1 ms Force decline time constant
    when strongly bound
    cross-bridges absent
    τ2 ms 124.4 Force decline time constant
    due to friction between actin
    and myosin
    Km Sensitivity of strongly bound
    cross-bridges to CN
    Arest Nms 3.009 Value of the variable A
    when muscle is not fatigued
    Km,rest 0.103 Value of the variable Km
    when muscle is not fatigued
    τ1,rest ms 50.95 The value of the variable τ1
    when muscle is not fatigued
    αA 1ms2 4.0107 Coefficient for the force-model
    variable A in the fatigue
    model
    αKm 1msN 1.9108 Coefficient for the force-model
    variable Km in the fatigue
    model
    ατ1 1N 2.1105 Coefficient for force-model
    variable τ1 in the fatigue
    model
    τfat s 127 Time constant controlling the
    recovery of (A,Km,τ1)

     | Show Table
    DownLoad: CSV

    For the purpose of our analysis the force-fatigue model is written as the single-input control system:

    dxdt(t)=F0(x(t))+u(t)F1(x(t)) (9)

    with x=(x1,x2,x3,x4,x5)=(CN,F,A,Km,τ1) where u denotes the control u=Es(t) corresponding to the sampled physical control data:

    (I1,I2,,In,η1,,ηn) (10)

    with constraints

    IminIiImax,0ηi1

    corresponding to interpulses bounds and amplitude convexification. This leads to a non-linear model with sampled control data with prescribed convex control constraints (η,I)C;η=(η1,,ηn),I=(I1,,In).

    The force model can be briefly investigated to describe preliminary results.

    Observe that in the force model the static non-linearities described by the mappings a and b can be absorbed by time reparameterization to provide an explicit form for the force responses.

    Proposition 1. The force model can be integrated by quadratures using a time reparameterization.

    Proof. Integrating the linear system (3), the equation (4) can be written as

    dFds(s)=c(s)F(s) (11)

    with ds=b(t)dt,c(s)=A(s)a(s)/b(s). It can be integrated using Lagrange formula.

    Thus this gives an explicit force response in the time parameter s

    (I,η,s)F(I,η,s). (12)

    Clearly we have

    Lemma 3.1. The above mapping is smooth with respect to I,η and piecewise smooth with respect to s.

    For the sake of providing a smooth response for the observer it is sufficient to smooth the physical sampled control date as follows: use a bump function to smooth each Heaviside mapping H(tti) at the sampling time ti.

    For the sake of the geometric analysis of the dynamics and to simplify the control constraints on (I,η) the FES signal Es(t) is taken as the input u() of the control system. Using (2) one can write

    u(t)=Es(t)=et/τcv(t)

    with

    v(t)=1τcni=1H(tti)Riηieti/τc. (13)

    Whence 0=t1<t2<<tn are given v(t) is a piecewise constant control depending upon the parameters η1,,ηn and the dynamics of the force model can be analyzed in the frame of geometric control.

    First of all the control system (9) is analyzed in the framework of geometric control.

    Consider a (smooth) control single-input control system.

    {dxdt=X(x)+uY(x)y=h(x) (14)

    where xRn,y=h(x)R corresponds to a (smooth) single-observation mapping. The following concepts rely on seminal results of geometric control, see [18,13].

    We denote [U,V] the Lie bracket of two (smooth) vector fields of Rn:

    [U,V](x)=Ux(x)V(x)Vx(x)U(x)

    and a vector field U acts on (smooth) mappings f with the Lie derivative:

    LUf=fxU(x).

    Let D1=span{X,Y} and define recursively : Dk=span{Dk[D,Dk1]},k>1. The length of [Xi1,,[Xik1,Xik]] is k. Hence Dk represents Lie brackets of X,Y with lengths smaller than k+1 and denote DL.A.=spank0Dk the Lie algebra generated by {X,Y}. The system is called weakly controllable if for each xRn,DL.A.(x)=Rn.

    The observation space is the set of mappings: Θ={LGh;GDL.A.} and the system is called observable if for each x1,x2Rn, x1x2 there exists θΘ such that θ(x1)θ(x2).

    Taking x0Rn, a frame at x0 is a set of elements X1,,Xn of DL.A. such that: X1,,Xn are linearly independent at x0 and ni=1length(Xi) is minimal.

    The system is called feedback linearizable in an open set VRn if ˙x=X(x)+uY(x) is feedback equivalent to the linear system ˙x=Ax+ub, that is there exists a diffeomorphism φ on V, φ(0)=0 and a feedback u=α(x)+β(x)v,β(x)0 such that g(X,Y)=(A,b) with g=(φ,α,β) acting by change of coordinates and (affine) feedback.

    Fix x(0)=x0 and T>0 and consider the extremity mapping: Ex0,T:u()L([0,T])x(T,x0,u) where x() is the response of ˙x=X(x)+uY(x) to u defined on [0,T].

    The control u() is called singular on [0,T] if the extremity mapping Ex0,T is not of maximal rank n when evaluated at u().

    Geometric analysis of an observed system of the form (14) amounts to compute DL.A., the observation space, the singular controls and the feedback equivalence properties. Achievements of geometric optimal control amounts to synthesize the optimal control in relation with the Lie brackets properties of a frame.

    The control in the force-fatigue model (1) falls into the framework of sampled-optimal control problem and we refer to [3] for more details. We use the following terminology.

    Consider a control system ˙x=f(x,u). When the state x() and the control u() evolve continuously in time, we speak of a continuous-time control problem and the control is said permanent. The sampled-data control case is when u() is piecewise constant.

    For simplicity consider the reduced model in dimension 3 related to the force-fatigue model.

    Let the minimization problem be defined by

    T0(u(s)2+(FFref)2)dsmin|u|1

    subject to the dynamics in ˜x=(CN,F,A) given by

    dCNdt=CNτc+udFdt=aAbFdAdt=AArestτfat+αAF (15)

    and with the initial conditions

    F(0)=A(0)=CN(0)=0.

    The parameters T,τ2,τfat,αA,Arest and the variables Km,τ1 are constant and fixed to some prescribed values.

    The system (15) can be written into the form

    ˙˜x=˜F0(˜x)+u˜F1(˜x).

    Remark 1. The model (15) is a simplification of the force-fatigue model. In particular, without loss of generality, we don't take into account the factor exp(t/τc)/τc appearing in (2). Control constraints |u|1 are not the physical constraints, see Section 3.3. Also the fatigue dynamics is reduced to a single equation, motivated by the controllability properties of the system (6)-(7)-(8).

    The permanent control case. The problem can be summarized as a permanent optimal control problem as follows

    {minu()T0(u(s)2+(FFref)2)ds˙˜x=˜F0(˜x)+u˜F1(˜x)u(t)[1,1](CN(0),F(0),A(0))=(0,0,0). (16)

    The pseudo-Hamiltonian of the system is

    H(˜x,p,p0,u)=p(˜F0(˜x)+u˜F1(˜x))+p0(u2+(FFref)2)

    where (p,p0):[0,T]R4 is the adjoint vector.

    We denote by Hi=p˜Fi,i=0,1 the Hamiltonian lifts of the vector fields ˜Fi,i=0,1.

    Normal case: p0=1/2. Applying the Pontryagin maximum principle (PMP), the optimal (permanent) control is given in the normal case by

    u(t)={1 if H1(˜x(t),p(t))1,1 if H1(˜x(t),p(t))1,H1(˜x(t),p(t))otherwise. (17)

    Abnormal case: p0=0. Abnormal controls are characterized by the equation p˜F1=0. Differentiating twice with respect to t leads to

    p[˜F0,˜F1]=0,p([[˜F0,˜F1],˜F0]+u[[˜F0,˜F1],˜F1])=0.

    Eliminating p we obtain D+uD=0 where

    D=det(˜F1,[˜F1,˜F0],[[˜F1,˜F0],˜F0])D=det(˜F1,[˜F1,˜F0],[[˜F1,˜F0],˜F1]).

    Computing we have

    D=αA(a(CN)Ab(CN)F)2,D=0

    and D=0 is equivalent to A(CN(τ1+τ2)+Kmτ1)2+Fτ2(CN+Km)2=0. This implies F=A=0 and the following lemma.

    Lemma 4.1. There are no physically admissible singular trajectories for the problem (15).

    The sampled-data control case. The corresponding optimal sampled-data control problem is given by

    {minT0(u(kTs)2+(FFref)2)dt,withk=t/Ts˙˜x=˜F0(˜x)+u˜F1(˜x)u(kTs)[1,1](CN(0),F(0),A(0))=(0,0,0). (18)

    where Ts>0 is a fixed sampling period such that T=jTs for some jN.

    Normal case: p0=1/2. Following [3], the optimal (sampled-data) control is

    u(kTs)argmaxy[1,1]1T(k+1)TskTsH(s,˜x(s),p(s),p0,y)ds (19)

    for all k{0,,d}.

    Write

    ˉH1=1T(k+1)TskTsp1(s)ds

    then, the optimal sampled-data control is

    u(kTs)={1 if ˉH1>1,1 if ˉH1<1,ˉH1 otherwise.  (20)

    Numerical results. In Fig. 1, we represent numerical results for several values of T/Ts. The optimal permanent control is represented with thin continuous line. A Gaussian quadrature rule is used to compute ˉH1. We observe the convergence of the optimal sampled-data control to the permanent control as Ts tends to 0.

    Figure 1. 

    Time evolution of the permanent control (thin continuous line) and sampled-data control for several values of the sampling period Ts{T/20,T/40,T/200}

    .

    The force evolution is compared for Km,rest and different values Km,rest (±30% of error) for I=10ms, 50ms and 100ms (see Fig. 2 for I=10ms). In the case of I=10ms, the maximum force error (see Fig. 3) is of 0.3%. Following interpulse value, the maximum force error is obtained for I=100ms (1.3%) which means that a tolerance of ±30% gives force evolution with a good accuracy.

    Figure 2. 

    Evolution of Km for different initial conditions (case of I=10ms)

    .
    Figure 3. 

    Relative error of the force for a well known and erroneous Km initial condition (case of I=10ms)

    .

    In this section, we design a modified version of the standard high-gain observer given in [11] taking into account the specific structure of the problem. The system defined by the force equation (4) and the fatigue model (6) and (8) can be rewritten as the single input-output system:

    {˙x(t)=am(t,Es(t))f(x(t),Es(t))y(t)=h(x(t))=F(t), (21)

    with x=(F,A,τ1)R3, yR, EsR, 0<a<1 is given by (5) and m is a positive integer. Note that in (21), Km is the solution of (7) thanks to the weak sensibility of the solution with respect to this variable (see Section 5.1 for the sensibility study of the force versus Km). We introduce the change of variables ϕ:

    {ϕ:R3R3xϕ(x)=[h(x),Lf(h(x)),L2f(h(x))]. (22)

    We have:

    ϕx=(1001/(τ1+τ2a)aF/(τ1+τ2a)2L2fh(x)x1L2fh(x)x2L2fh(x)x3)

    where

    L2fh(x)x1=aαA+1(x3+τ2a)2(1x3τ1,restτfat+2ατ1x1),L2fh(x)x2=aτfat+x3+τ2aτfat(x3+τ2a),L2fh(x)x3=1(x3+τ2a)2(x2ax1τfat)2x1(x3+τ2a)3(1x3τ1,restτfat+ατ1x1)

    and

    det(ϕx)=a(x2aτfat(x3+τ2a)τfatx1+2x1(x3τ1,restατ1τfatx1)τfat(x3+τ2a)3 (23)

    A sensibility study in (6) and (8) concerning respectively (AArest)/τfat and (τ1τ1,rest)/τfat, shows that neglecting these two terms induces a maximum error of 7% for A and τ1, and 6% for the force value. Hence we use the simplified model:

    dAdt=αAF (24)
    dτ1dt=ατ1F. (25)

    Hence (4) gives:

    det(ϕx)=a(x2a(x3+τ2a)x12ατ1x21)(x3+τ2a)3. (26)

    Since Es(t)0, we have:

    t[0,T]:(CN,x1)0,(x2,x3,Km)>0(x3+τ2a)3>0.

    The equation det(ϕx)=0 yields

    a=0CN=0(CN=0x1=0)(period of rest) (27)
    x2a(x3+τ2a)x12ατ1x21=0x1=11+8ατ1x2a(x3+τ2a)4ατ1. (28)

    Since ατ1105, 0a1 and (x3+τ2a)102, we have from (28)

    8ατ1x2a(x3+τ2a)1(x1=0(CN=0),x1=12ατ1).

    Thus (ϕx(ˆx(t))=0,ˆx=(ˆF,ˆA,ˆτ1) for:

    (CN,x1)=0 which corresponds to the period of rest,

    x1=1/(2ατ1), x1105 which is greater than the maximum force value.

    In the simplified model (4), (24) and (25), ϕx(ˆx(t)) is invertible during stimulation period (CN,x1>0), and the matrix is well-conditioned.

    Based on the simplified model, the modified high-gain observer is defined as:

    ˙ˆx(t)=amf1(ˆx(t),Es(t))am(ϕx(ˆx(t))1S1θCT(Cˆx(t)y(t)). (29)

    The experimental choice of m depends on the pulses frequency. Sθ is a symmetric definite positive matrix given by the following Lyapunov equation:

    θSθ(t)+ATSθ(t)+Sθ(t)A=CTC (30)

    where θ is a high-gain tuning parameter introduced in [11], and

    A=(010001000),C=(100).

    Hence the components of Sθ=[Sθ(l,k)]1l,k3 have the following form:

    Sθ(l,k)=(1)l+k(l+k2k1)θ(l+k1),(nk)=n!(nk)!k!. (31)

    Let

    ˙x(t)=am(t)f1(x(t),Es(t))=f(x(t),Es(t))y(t)=h(x(t)) (32)

    and

    z1=h(x),˙z1=Lfh(x)=amLf1h(x)=amz2,˙z2=Lf(Lf1h(x))=amLf1(Lf1h(x))=amz3,˙zn1=amzn,˙zn=amLnf1h(x)=φn(u,z).

    Then

    (˙z1˙z2˙zn)=am(0100100)(z1z2zn)+(00φn(u,z)). (33)

    Therefore

    {˙z=amAz+φ(u,z)y=Cz. (34)

    Take C=(100) and consider the system (33). We make the following assumptions.

    ● H1) φn(u,z) is globally Lipschitz with respect to z and uniformly with respect to u,

    ● H2) Let U be the control domain, a compact KRn and two positive constants amin and amax such that: for all u valued in U and initial conditions z(0) K, we have amina(t)amax,

    ● H3) a() is C1 (see Remark 2 below),

    ● H4) m˙a(t)/a(t)<1.

    The observer has the following expression:

    ˙ˆz(t)=a(t)mAˆz(t)a(t)mS1θCT(Cˆz(t)y(t)) (35)

    where Sθ is the solution of the Lyapunov equation:

    θSθ+ATSθ+SθACTC=0. (36)

    The solution of (36) can be rewritten as

    Sθ=1θDθS1Dθ. (37)

    and

    Dθ=diag(1,1θ,,1θn) (38)

    with S1 is the solution of (36) for θ=1. Let ρ=am and e=ˆzz ˙e=˙ˆz˙z, then

    ˙e=ρAˆz+ˆφ(u,ˆz)S1θCT(Cˆz(t)y(t))ρAzφ(u,z)=ρ(AS1θCTC)e+ˆφφ (39)

    Write ˉe=ρDθe, then

    ˙ˉe=˙ρDθe+ρDθ˙e=˙ρDθe+ρDθ(ρ(AS1θCTC)e+ˆφφ)=˙ρDθe+ρDθ(ρ(AθD1θS11D1θDθCTCDθ)e+ˆφφ),DθCTCDθ=CTC=˙ρρˉe+ρθAˉeθS11CTCˉe+ρDθ(ˆφφ),DθAD1θ=θA (40)

    Consider the Lyapunov function: V=ˉeTS1ˉe. Then

    ˙V=2ˉeTS1˙ˉe=2ˉeTS1(˙ρρˉe+ρθAˉeθS11CTCˉe+ρDθ(ˆφφ)) (41)

    2ˉeTS1Aˉe=ˉeT(S1A+ATS1)ˉe and S1A+ATS1=S1+CTC.

    Then, we have

    ˙V=θρV+θρˉeTCTCˉe2ˉeTCTCˉe+2ˉeTS1˙ρρˉe+2ˉeTS1ρDθ(ˆφφ)=θρV+(θρ2)∣∣Cˉe2+2˙ρρV+2ˉeTS1ρDθ(ˆφφ)=(θρ2˙ρρ)V+(θρ2)∣∣Cˉe2+2ˉeTS1ρDθ(ˆφφ) (42)

    We deduce that

    θρ2˙ρρ>0θ>2˙ρρ2,θρ2<0θ<2ρ2˙ρρ2<θ<2ρ
    (˙ρρ=m˙aa)<1. (43)

    Using (43) and assumption H1:

    ˙VγV+2ˉeTS1ρDθ∣∣ˆφφ∣∣,γ>0γV+2ˉeTS1ρθn^φnφnγV+2ˉeTS1ρθnK1ˉe(γ2amθnK1)V. (44)

    For m and θ sufficiently large: γ>2am/θnK1˙Vγ1V,γ1>0.

    Remark 2. In the force-fatigue model, a(t) is piecewise smooth, the lack of regularity is numerically bypassed by the choice of the integer m. For example, for I=10ms, m=3 is sufficient to estimate the whole variables. However, for I=25ms, m must be at least equal to 5 (see observer simulations in Section 7).

    MPC computes a sequence of predicted controls to optimize plant behavior (see Fig. 4).

    Figure 4. 

    General MPC strategy diagram

    .

    Model Predictive Heuristic Control (MPHC) was introduced by Richalet et al. [17] using impulse response type dynamic model. Dynamic Matrix Control (DMC) followed in 1980 (Cutler et al. [5]) using step response type dynamic model. State space formulation of MPC was introduced by Li et al. [16].

    Consider a system written as

    ˙x=f(x,u),x(0)=x0 (45)

    and a general cost function to be minimized

    J(u(.),x0)=0q(xu(τ,x0),u(τ))dτ (46)

    with

    xu(t,x0)=x0+t0f(xu(τ,x0),u(τ))dτ (47)

    and

    q(0,0)=0,q()C2q(x,u)cq(||x||2+||u||2),cq>0uq(x,u) is convex for all x. (48)

    Without constraints, Bellman's principle of optimality (1950) gives the solution. In the case of constrained problem, instead of J, we minimize the following cost:

    J(u(.),x0,T)=T0q(xu(τ,x0),u(τ))dτ+v(xu(T,x0)) (49)

    with

    J(x0,T)=minu(.)J(u(.),x0,T) (50)

    and

    u(t,x0,T)=argminu(.)J(u(.),x0,T). (51)

    This optimization over a finite horizon follows the algorithm:

    1. Solve minJ(u(.),x0,T) and find u(.,x0,T).

    2. Apply u(.,x0,T) for τ[0,Ts[, 0TsT, (Ts: sampling period).

    3. Repeat using x(Ts) instead of x0.

    To explain the method, we consider the discrete linear system:

    {x(k+1)=Ax(k)+Bu(k)y(k)=Cx(k) (52)

    where x(t) Rn, u(t) Rm, y(t) Rp are respectively the state, input and output vectors, and t=kTs. This system (66) results from a discrete system modelling or continuous-discrete transformation. Practical considerations of implementation make this form more appropriate in the framework of MPC.

    We suppose that the system is both controllable and observable, and we denote by x(k+i/k),i0, the predicted state vector x(k+i) with initial condition x(kTs).

    Consider the unconstrained problem and the quadratic cost:

    J(k)=j=k[(y(j+1/k)yref)TQ(y(j+1/k)yref)+(u(j/k)uref)TR(u(j/k)uref)+Δu(j/k)TSΔu(j/k)] (53)

    with Q=QT 0, R=RT 0 and S=ST 0. The couple (yref, uref) is solution of (52) and Δu(k)=u(k)u(k1). The solution of this problem is obtained using the LQ controller.

    In the case of constrained problem, LQR could not be used to solve control problem. System (52) allows to calculate x(k+i/k), i=1,...,Nr at time kTs, (Nr being the receding horizon length):

    x(k+i/k)=Aix(k/k)+i1j=0Aij1Bu(k+j/k). (54)

    The input-output relation along the receding horizon is then:

    (y(k+1/k)y(k+2/k)y(k+Nr/k))=(CACA2CANr)x(k/k)+(CB00CABCB0CANr1BCABCB)(u(k/k)u(k+1/k)u(k+Nr1/k)). (55)

    Write (55) as

    ˉy=Ψx(k/k)+Γˉu (56)

    where ˉyk RpNr, ˉuk RmNr, Ψ RpNr×n, Γ RpNr×mNr.

    Constraints can take the form:

    ● Control variable:

    uminu(k)umax. (57)

    ● Change rate of the control variables:

    ΔuminΔu(k)Δumax (58)

    with

    Δu(k)=u(k)u(k1). (59)

    ● Soft output variables constraints (relaxed output constraints using large slack variable sv to avoid constraints conflicts when solving control problem):

    yminsv1y(k)ymax+sv1. (60)

    ● Soft state variables constraints (for the same reason as output constraints):

    xminsv2x(k)xmax+sv2 (61)

    The cost function (53) becomes:

    J(k)=k+Nr1j=k[(y(j+1/k)yref)TQ(y(j+1/k)yref)+(u(j/k)uref)TR(u(j/k)uref)+Δu(j/k)TSΔu(j/k)] (62)

    subject to (56) and a set of constraints among (57)-(61).

    This minimization problem takes the form:

    min1/2XTEX+XTFsubject to ˉy=Ψx(k/k)+ΓˉuMXγ (63)

    and ˉy=Cˉx with X=[ˉx,ˉu]T. Optimization problem (63) is a convex Quadratic Programming (QP) problem. Denote X=[ˉx,ˉu]T the global minimizer at each iteration. Once ˉu is calculated, only u(k/k) is applied at time t=kTs and we iterate the algorithm.

    In the case of NMPC (non-linear Model Predictive Controller), the dynamics is defined using a discretization of the equation (45). Using (62), optimization problem becomes:

    min1/2XTEX+XTFsubject to x(k+i/k)=g(x(k+j/k),u(k+j/k)),j=0,,Nr1MXγ (64)

    with X=[ˉx,ˉu]T. Equality constraint in the optimization problem (64) is non-linear. Now convexity condition in (64) is not guaranteed, then X could be a local minimizer. NMPC (non-linear Model Predictive controller) is used to solve this problem.

    Numerical solutions are computed using Active Set, Primal-Dual or SQP methods (Wang [19], Fletcher [10] and Boyd and Vandenberghe [4]).

    In the case of force-fatigue model, one must consider the more general case of a varying Ts (Ts I(i)), then u(i)=[I(i)η(i)]T, where I(i)=t(i)t(i1) is the interpulse between two successive pulses and η(i) is the pulse amplitude applied at time ti. In this case, ˉu(k)=[u(k/k)u((k+Nr1)/k)] is the discrete receding horizon control vector (of length Nr) to be calculated by the NMPC at time t corresponding to iteration k. Then:

    ˉu(k)=(I(k/k)I(k+1/k)I((t(Nr))/k)η(tk/k)η((tk+I(k/k))/k)η((tk+Nr1j=0I(k+j))/k)) (65)

    t(Nr) being the final time of the optimization horizon which is unknown a priori. For instance, using single move strategy, we have u(k/k)=u(k+1/k)==u((k+Nr1)/k), (see Fig. 5).

    Figure 5. 

    (left half plane) Es and force profile for applied amplitude and interpulse stimulation, (right half plane) Predicted Es and force using single move strategy to be optimized

    .

    The force-fatigue model being non-linear, an NMPC is used to solve the problem with Non-linear Programming Algorithm (NPA). In this case, criterion to be minimized is:

    J(k)=k+Nr1j=k(F(j+1/k)Fref)2 (66)

    subject to: 0η(i)1 and 0.01msI(i)0.1ms.

    The estimation of the state variables vector (by the high-gain observer) is used as an initial variables vector to perform the NMPC over the horizon Nr.

    To exhibit the interest of the use of the power m in the non-linear observer, consider a MPC based on a non-linear observer with only amplitude as control variable. We consider also the worst case of +30% of error of Km.

    The following simulation results follow the protocol used in practice (a set of periods with stimulation and rest time slots in each period). During the stimulation time slot, the control is calculated to bring the force to Fref and in the rest time slot, the stimulation amplitude is set to 0. In this section, only two stimulation periods are considered.

    For I=10ms (100Hz) and I=25ms (40Hz), we take m=3 and m=5 respectively. Fig. 6 and Fig. 7 represent the ˆA for I=10ms and I=25ms, respectively. Fig. 8 and Fig. 9 are the ˆτ1 for I=10ms and I=25ms. ˆA converges after 50ms when I=10ms and 100ms when I=25ms. Concerning ^τ1, it converges after 75ms when I=10ms and 200ms when I=25ms. Large I seems to delay the convergence of the observer.

    Figure 6. 

    Evolution of A and ˆA for I=10, 30% error of Km

    .
    Figure 7. 

    Evolution of A and ˆA for I=25, 30% error of Km

    .
    Figure 8. 

    Evolution of τ1 and ^τ1 for I=10, 30% error of Km

    .
    Figure 9. 

    Evolution of τ1 and ^τ1 for I=25, 30% error of Km

    .

    Fig. 10 represents the force response for amplitude control strategy (for a receding horizon Nr=10) based on the proposed observer for I=25ms and a force reference of 250N. Force mean value converges to the force reference after 200ms.

    Figure 10. 

    Evolution of F, ˆF and F mean value over I for I=25, 30% error of Km, Fref=250N

    .

    In this section, five stimulations periods are presented (due to the experimental protocol). Figures 11, 12, and 13, are the force response in the case of Fref=425N and (Nr=3,5,10), the interpulse and the amplitude controls for Nr=10, respectively. As expected, Nr=10 gives the best regulation performances (response time and overshoot).

    Figure 11. 

    Evolution of the force for a reference force of 425N and different receding horizons (3,5 and 10)

    .
    Figure 12. 

    Evolution of the interpulse (control) for a reference force of 425N and a preditive horizon of 10

    .
    Figure 13. 

    Evolution of the amplitude (control) for a reference force of 425N and a preditive horizon of 10

    .

    Before t=6000ms with Nr=10, the force is correctly maintained at Fref=425N. For t6000ms (starting from the forth period), the fatigue level A is very high so that the maximum value of the amplitude control (see Fig. 13) and the interpulse control (see Fig. 12) cannot maintain F at Fref. Maximum interpulse frequency is not used at the beginning of the forth and the fifth periods. Higher value of Nr could correct this problem (supposing that global minimizer of MPC algorithm is reached at each iteration). However, increasing Nr will render MPC algorithm time consuming and cause a problem for real time implementation.

    Figure 14. 

    Evolution of the interpulse (control) for a reference force of 425N and a preditive horizon of 10

    .
    Figure 15. 

    Evolution of the amplitude (control) for a reference force of 425N and a preditive horizon of 3

    .

    This work deals with the control and estimation of the state variables of the Ding et al. force-fatigue model where the control is the interpulse and/or amplitude of the electrical stimulation. Preliminary geometric analysis of the force control controlling force level is presented with the aim of a future PMP control strategy. In the case of a fixed interpulse, the proposed high-gain observer using the force measurements exhibits the relation between the interpulse and the parameter m to perform accurate variables estimation. Model Predictive Control (MPC) strategy is presented in the case of both stimulation interpulse and amplitude as control variables. Simulation results are presented concerning the high-gain observer and the MPC force control. These simulations show the effect of the receding horizon on the control efficiency. Reasonable value of Nr is however suitable to guarantee a short computation time for real time implementation.



    [1] J. L. Lions, On some questions in boundary value problems of mathematical physics, North-Holland: North-Holland Math. Stud., 30 (1978), 284–346. https://doi.org/10.1016/S0304-0208(08)70870-3 doi: 10.1016/S0304-0208(08)70870-3
    [2] H. Chen, M. M. Fall, B. Zhang, On isolated singularities of Kirchhoff equations, Adv. Nonlinear Anal., 10 (2021), 102–120. https://doi.org/10.1515/anona-2020-0103 doi: 10.1515/anona-2020-0103
    [3] W. He, D. Qin, Q. Wu, Existence, multiplicity and nonexistence results for Kirchhoff type equations, Adv. Nonlinear Anal., 10 (2021), 616–635. https://doi.org/10.1515/anona-2020-0154 doi: 10.1515/anona-2020-0154
    [4] A. Hamydy, M. Massar, N. Tsouli, Existence of solutions for p-Kirchhoff type problems with critical exponent, Electron. J. Differ. Equ., 105 (2011), 1–8.
    [5] E. Dibenedetto, Degenerate Parabolic Equations, Springer, New York, 1993. https://doi.org/10.1007/978-1-4612-0895-2
    [6] H. Ishii, Asymptotic stability and blowing up of solutions of some nonlinear equations, J. Differ. Equ., 26 (1977), 291–319. https://doi.org/10.1016/0022-0396(77)90196-6 doi: 10.1016/0022-0396(77)90196-6
    [7] M. Tsutsumi, Existence and nonexistence of global solutions for nonlinear parabolic equations, Publ. Res. Inst. Math. Sci., 8 (1972), 211–229. https://doi.org/10.2977/prims/1195193108 doi: 10.2977/prims/1195193108
    [8] H. A. Levine, L. E. Payne, Nonexistence of global weak solutions of classes of nonlinear wave and parabolic equations, J. Math. Anal. Appl., 55 (1976), 329–334. https://doi.org/10.1016/0022-247X(76)90163-3 doi: 10.1016/0022-247X(76)90163-3
    [9] M. Ghisi, M. Gobbino, Hyperbolic-parabolic singular perturbation for mildly degenerate Kirchhoff equations: time-decay estimates, J. Differ. Equ., 245 (2008), 2979–3007. https://doi.org/10.1016/j.jde.2008.04.017 doi: 10.1016/j.jde.2008.04.017
    [10] Q. Lin, X. Tian, R. Xu, M. Zhang, Blow up and blow up time for degenerate Kirchhoff-type wave problems involving the fractional Laplacian with arbitrary positive initial energy, Discrete Contin. Dyn. Syst. Ser. S, 13 (2020), 2095–2107. https://doi.org/10.3934/dcdss.2020160 doi: 10.3934/dcdss.2020160
    [11] N. Pan, P. Pucci, R. Xu, B. Zhang, Degenerate Kirchhoff-type wave problems involving the fractional Laplacian with nonlinear damping and source terms, J. Evol. Equ., 19 (2019), 615–643. https://doi.org/10.1007/s00028-019-00489-6 doi: 10.1007/s00028-019-00489-6
    [12] X. Wang, Y. Chen, Y. Yang, J. Li, Kirchhoff-type system with linear weak damping and logarithmic nonlinearities, Nonlinear Anal., 188 (2019), 475–499. https://doi.org/10.1016/j.na.2019.06.019 doi: 10.1016/j.na.2019.06.019
    [13] M. Chipot, T. Savitska, Nonlocal p-Laplace equations depending on the Lp norm of the Gradient, Adv. Differ. Equ., 19 (2014), 997–1020.
    [14] Y. Han, Q. Li, Threshold results for the existence of global and blow-up solutions to Kirchhoff equations with arbitrary initial energy, Comput. Math. Appl., 75 (2018), 3283–3297. https://doi.org/10.1016/j.camwa.2018.01.047 doi: 10.1016/j.camwa.2018.01.047
    [15] S. Zheng, M. Chipot, Asymptotic behavior of solutions to nonlinear parabolic equations with nonlocal terms, Asymptotic Anal., 45 (2005), 301–312.
    [16] Y. Fu, M. Xiang, Existence of solutions for parabolic equations of Kirchhoff type involving variable exponent, Appl. Anal., 95 (2016), 524–544. https://doi.org/10.1080/00036811.2015.1022153 doi: 10.1080/00036811.2015.1022153
    [17] J. Li, Y. Han, Global existence and finite time blow-up of solutions to a nonlocal p-Laplace equation, Math. Model. Anal., 24 (2019), 195–217. https://doi.org/10.3846/mma.2019.014 doi: 10.3846/mma.2019.014
    [18] L. E. Payne, D. H. Sattinger, Saddle points and instability of nonlinear hyperbolic equtions, Israel J. Math., 22 (1975), 273–303. https://doi.org/10.1007/BF02761595 doi: 10.1007/BF02761595
    [19] R. Xu, J. Su, Global existence and finite time blow-up for a class of semilinear pseudo-parabolic equations, J. Funct. Anal., 264 (2013) 2732–2763. https://doi.org/10.1016/j.jfa.2013.03.010
    [20] H. A. Levine, Some nonexistence and stability theorems for solutions of formally parabolic equations of the form Put=Au+F(u), Arch. Ration. Mech. Anal., 51 (1973), 371–386. https://doi.org/10.1007/BF00263041 doi: 10.1007/BF00263041
    [21] M. Liao, W. Gao, Blow-up phenomena for a nonlocal p-Laplace equation with Neumann boundary conditions, Arch. Math., 108 (2017), 313–324. https://doi.org/10.1007/s00013-016-0986-z doi: 10.1007/s00013-016-0986-z
    [22] H. Brezis, Functional Analysis, Sobolev spaces and partial differential equations, Springer, New York, 2010.
    [23] G. A. Philippin, V. Proytcheva, Some remarks on the asymptotic behaviour of the solutions of a class of parabolic problems, Math. Methods Appl. Sci., 29 (2006), 297–307. https://doi.org/10.1002/mma.679 doi: 10.1002/mma.679
    [24] Y. Han, Finite time blowup for a semilinear pseudo-parabolic equation with general nonlinearity, Appl. Math. Lett., 99 (2020), 1–7. https://doi.org/10.1016/j.aml.2019.07.017 doi: 10.1016/j.aml.2019.07.017
    [25] Y. Li, C. Xie, Blow-up for p-Laplacian parabolic equations, Electron. J. Differ. Equ., 20 (2003), 1–12.
  • This article has been cited by:

    1. Toufik Bakir, Bernard Bonnard, Loïc Bourdin, Jérémy Rouot, Pontryagin-Type Conditions for Optimal Muscular Force Response to Functional Electrical Stimulations, 2020, 184, 0022-3239, 581, 10.1007/s10957-019-01599-4
    2. Bernard Bonnard, Jérémy Rouot, Geometric optimal techniques to control the muscular force response to functional electrical stimulation using a non-isometric force-fatigue model, 2021, 13, 1941-4897, 1, 10.3934/jgm.2020032
    3. Toufik Bakir, Bernard Bonnard, Sandrine Gayrard, Jérémy Rouot, Finite dimensional approximation to muscular response in force-fatigue dynamics using functional electrical stimulation, 2022, 144, 00051098, 110464, 10.1016/j.automatica.2022.110464
  • Reader Comments
  • © 2022 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(1818) PDF downloads(120) Cited by(0)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog