Tutorial #4: Seeking optimal strokes

Biohydrodynamics Toolbox    

Tutorial #4: Seeking optimal strokes

In this page we explain through an example how a swimming strategy (actually in our example it will rather be here a flying strategy) can be easily optimized using the built-in MATLAB function fminsearch.

What means getting a stroke optimized?

In the absence of fluid (for example by setting fluid density = 0 in the DAT-File) there would be of course no hydrodynamical forces applied to the bodies. Since the deformations are due to inner forces only, the articulated bodies would not be able to move their centers of mass neither to modify their angular momenta. It is a straight consequence of Newton's first and second laws. Thus, the bodies' locomotion results from the interaction with the surrounding fluid. The deformations of the bodies set the fluid into motion which in return (Newton's third law) pushes the bodies. The locomotion occurs due to the transfer of momentum between the bodies and the fluid. Now the question we are concerned with is: how to quantify the efficiency of a locomotion strategy? We dispose for that of three physical amounts: energy, power and torque.

Energy

The fluid model BhT deals with is frictionless ; there is no energy dissipation. Nevertheless, by running the function bht_energy, one can observe that sometimes fluid-bodies systems undergo variations of the overall energy. These variations are only due to what happen inside the bodies (i.e. the bodies' inner Dynamics). Indeed, imagine for a while that all of the immersed bodies' hinges behave like springs and that the potential energy of these springs be taken into account in the energy balance of the system. Then, the locomotion strategy would be energetically neutral.  The same conclusion can be drawn replacing the springs by engines that would alternatively be able to work as generators and then to store back energy. For these cases, energy is definitively not a relevant quantity to measure the efficiency of the locomotion strategy.

However, as far as common engines are concerned (which are not able to store back energy), they expend energy as soon as the system fluid-bodies undergo variations of the overall energy. The physical amount that measures energy variations is the power.

Power

The power expended by the bodies at its joints is obtained by multiplying the torques by the angular velocities. Since the fluid model is frictionless, the total power expended by all of the joints of all of the immersed bodies coincides with the time derivative of the total energy of the system fluid-bodies (this can be checked easily using the function bht_energy). In other words, any variation of the overall energy results from what happen at the joints. 

A classical notion of efficiency for fish locomotion, based on powers considerations, is due to Lighthill (1952). It is defined by the inverse of the ratio between the average power expended by the body and the power that an external force would spend to produce the same displacement. Unfortunately, this notion is meaningless in a perfect fluid since the average power expended by a bodies over a stroke can be zero. Indeed, let us consider an articulated body alone in a fluid of infinite extend, performing periodic strokes. The total energy of the system fluid-body is also periodic, with the same period as the strokes and hence the average power over one stroke is null. This is due to the fact that all of the positions the body occupies in the fluid are equivalent (because the fluid is of infinite extend) and hence that the total energy depends only on its shape and velocity (this is no longer true in the presence of several bodies or fixed boundaries). In fact, we are faced again with the same paradox described in the preceeding section but expressed here in terms of power rather than energy.

As we pointed it out in the preceeding section, such a reasoning presupposed the bodies to be able to store back energy in its course and it is not the case for common engines neither for muscles.  Therefore, one can rather try to get the stroke optimized over the time interval [Tmin Tmax] by minimizing the quantity:

 

where we use the notation introduced in the page bht_energy. 

Torque

The following example proves that optimizing only in terms of power has limitation especially when articulated bodies model animals. Consider a body composed of three ellipse-shaped solids connected together in line. If the hinges are locked (both angles = 0 for all time), the body behaves in the fluid like a rigid body. We launch it at the time t = 0 with non zero initial velocity and non zero initial angular velocity (click on the image to play the animation):

We next study how the torques and the power expended by the joints evolve in time by running bht_energy. We obtain the following results:

The powers of the joints remain always zero (because the joints have always null angular velocities) but the torques don't. If the articulated body models an artificial mechanism, we can assume that it is able to lock its joints with zero energy cost. However, if the articulated body models an animal, its muscles do work during the course ! In this case, getting the stroke optimized leads to minimize the amount:

Let us now treat an example thoroughly.

Experiment's features

The experiment consists of a flapping neutrally buoyant bird-like articulated body in an unbounded fluid. The bird is composed of a body with two wings, each one composed of three ellipse-shaped solids. The DAT-File is as follows (you can copy-past it into the MATLAB editor and save it as 'bird.dat').

fluid density = 1 mesh size = 0.1 max time= 6*pi time step = 0.1 fish =  {     initial position = [0;0;0]     initial velocity = [0;0;0]     colormap = winter   link = { label = body            mfilename = bht_ellipse            settings = [0 0 1 1 1] }    link = { label = rwing1            father = body            mfilename = bht_ellipse            settings = [0 0 1.5 0.3 1]            hinge coord = [1.2;0]            hinge local coord = [-1.7,0] }   link = { label = rwing2            father = rwing1            mfilename = bht_ellipse            settings = [0 0 1.5 0.3 1]            hinge coord = [1.7;0]            hinge local coord = [-1.7,0] }   link = { label = rwing3            father = rwing2            mfilename = bht_ellipse            settings = [0 0 1.2 0.3 1]            hinge coord = [1.7;0]            hinge local coord = [-1.4,0] }   link = { label = lwing1            father = body            mfilename = bht_ellipse            settings = [0 0 1.5 0.3 1]            hinge coord = [-1.2;0]            hinge local coord = [1.7,0]  }  link = { label = lwing2           father = lwing1           mfilename = bht_ellipse           settings = [0 0 1.5 0.3 1]           hinge coord = [-1.7;0]           hinge local coord = [1.7,0]   }  link = { label = lwing3           father = lwing2           mfilename = bht_ellipse           settings = [0 0 1.2 0.3 1]           hinge coord = [-1.7;0]           hinge local coord = [1.4,0]   } }

The controls M-File

A controls M-File is used to drive the six hinges' angles.  You can generate a pattern for this file by entering the command line:

>> bht_data_compile('DataFilename','bird','ControlsFilename','bird_cont');

The process of compiling the DAT-File must exit free of errors (but with several Warnings corresponding to the omitted optional fields). The pattern of the controls M-File bird_cont need to be next completed/modified as follows:

function [c,dc,dttc] = bird_cont(t,cont_parameters) % ====================================================== % preallocating output variables c0 = zeros(6,1); dc0 = zeros(6,1); dttc0 = zeros(6,1); % ====================================================== a = cont_parameters(1:3); ampl = cont_parameters(4:6); alpha(1:2) = cont_parameters(7:8); % ====================================================== % Warning: non zero controls' velocities (dc) at the % time t = 0 may produce fish's strong drift. % ====================================================== % % % ====================================================== %                     fish 1 % ====================================================== %              hinge              control number % ______________________________________________________ %          body-rwing1                1 %        rwing1-rwing2                2 %        rwing2-rwing3                3 %          body-lwing1                4 %        lwing1-lwing2                5 %        lwing2-lwing3                6 % ====================================================== [y,dty,dtty] = beat([0,pi/2,6*pi-pi/2,6*pi,6*pi],t); % ====================================================== c0(1) = ampl(1)*sin(t-a(1)); dc0(1) = ampl(1)*cos(t-a(1)); dttc0(1) = -ampl(1)*sin(t-a(1)); % ====================================================== c0(2) = ampl(2)*(sin(t-a(2))-alpha(1)); dc0(2) = ampl(2)*cos(t-a(2)); dttc0(2) = -ampl(2)*sin(t-a(2)); % ====================================================== c0(3) = ampl(3)*(sin(t-a(3))-alpha(2)); dc0(3) = ampl(3)*cos(t-a(3)); dttc0(3) = -ampl(3)*sin(t-a(3)); % ====================================================== c0(4) = -c0(1); dc0(4) = -dc0(1); dttc0(4) = -dttc0(1); % ====================================================== c0(5) = -c0(2); dc0(5) = -dc0(2); dttc0(5) = -dttc0(2); % ====================================================== c0(6) = -c0(3); dc0(6) = -dc0(3); dttc0(6) = -dttc0(3); % ====================================================== c = y*c0; dc = (y*dc0+dty*c0); dttc = (dttc0*y+2*dty*dc0+c0*dtty); % ====================================================== % Use the fonction BHT_CONTROLS_CHECK to check the formula

The flapping regime, from rest to rest, can be decomposed into three phases: the starting regime, from time t = 0 up to t = pi/2, the cruising regime from t = pi/2 up to t = 6*pi-pi/2 and the stopping regime from t = 6*pi-pi/2 up to t = 6*pi. To split the regime into three phases, we first set up the auxiliary variables c0, dc0 and dttc0 and next multiply them by the values returned by a function named beat to eventually get the output variables c, dc, dttc. The function beat returns the values of a smooth (C³), periodic function together with its first and second derivatives computed at the point t. The function beat has an extra input variable alpha which is a parameters array of five elements, alpha(5) being the period  (see graph below for the meaning of the other parameters). 

Starting and stopping regimes are introduced because the bird-like body has to start from rest to avoid a unwanted drift in its motion. Here is the script of the function beat you can copy-past into the MATLAB editor. Save it as beat.m.

function [x,dx,dttx] = beat(alpha,t) %------------------ t = mod(t,alpha(5)); %------------------------------------------------------ P = [35/128,0,-45/32,0,189/64,0,-105/32,0,315/128,0]; Pp = [315/128,0,-315/32,0,945/64,0,-315/32,0,315/128]; Ppp = [315/16,0,-945/16,0,945/16,0,-315/16,0]; %------------------------------------------------------ if t <= alpha(1);     x = 0;     dx = 0;     dttx = 0; elseif t<= alpha(2)     t1 = 2*(t-(alpha(2)+alpha(1))/2)./(alpha(2)-alpha(1));     %---------------------------------------------------     x = 0.5*polyval(P,t1)+0.5;     dx = polyval(Pp,t1)./(alpha(2)-alpha(1));     dttx = 2*polyval(Ppp,t1)./(alpha(2)-alpha(1))^2;     %--------------------------------------------------- elseif t<= alpha(3)     x = 1;     dx = 0;     dttx = 0;     %--------------------------------------------------- elseif t<= alpha(4);     t1 = 2*(t-(alpha(3)+alpha(4))/2)./(alpha(4)-alpha(3));     x = - 0.5*polyval(P,t1)+0.5;     dx = - polyval(Pp,t1)./(alpha(4)-alpha(3));     dttx = - 2*polyval(Ppp,t1)./(alpha(4)-alpha(3))^2;     %--------------------------------------------------- else     x =0;     dx =0;     dttx =0; end

The input variable cont_parameters in the controls M-File bird_cont is a vector of height elements entering in the definitions of the angles' controls c(k). You can perform some tests by running the function bht_kine_check and observe how modifying these parameters impact the bird's flapping. Try for instance:

>> bht_kine_check('DataFilename','bird','ControlsFilename','bird_cont','ControlsParameters',[1 2 3 pi/4 pi/6 pi/6 1 1]);

Here are some screenshots of what you might obtain (click on the first picture to see the animation):

Remark that whatever the variable cont_parameters is, the deformations are always 2*pi periodic within the cruising regime.

Computation of the bird's motion

Enter now the command line:

>> bht_traject_compute('datafilename','bird','controlsfilename','bird_cont','controlsparameters',[1 2 3 pi/4 pi/6 pi/6 1 1]);

After a while, you get the solution saved in the MAT-File bird_results. You can have a look at what the flapping looks like by using the function bht_simulation with the options 'axis', 'displaytime' and 'drawfishaxis' as follows:

>> film = bht_simulation('datafilename','bird','axis',[-11 11 -5 13],'displaytime','on','drawfishaxis',{3,'-r'});

Observe that although the bird's body undergo oscillations, its center of mass goes up in a much more regular way. Click on the image to play the movie: 

First optimization process

We wish now to adjust the variable cont_parameters to get the bird to go higher. To achieve this task, we use the following function altitude:

function height = altitude(cont_parameters) % ========================================================================= ampl = cont_parameters(4:6); alpha(1:2) = cont_parameters(7:8); % test overlapping ======================================================== recov = 1.2+3.4*cos(ampl(1))+3.4*cos(ampl(1)+ampl(2)*(1+alpha(1)))...     +2.6*cos(ampl(1)+ampl(2)*(1+alpha(1))+ampl(3)*(1+alpha(2))); % computations ============================================================ if recov > 0.1     [T,p,dp,q,dq,Q] = bht_traject_compute('datafilename','bird',...         'controlsfilename','bird_cont',...         'controlsparameters',cont_parameters,...         'saveresults','off');     %------------------     height = -Q(2,end);     %------------------ else     %----------     height = 0;     %---------- end

To avoid solids' overlaps, we compute the position (= recov)  of the right wing's extremity when the wings are down. If recov > 0, the wings will not overlap in their courses (we test rather recov > 0.1 for more security). If cont_parameters makes the wings overlap, the function returns 0. If not, the function returns the opposite of the maximum height reached by the bird since we have in mind to use the function fminsearch which only seeks minima of functions.

Thus, we run then the script:

% initial guess param_init = [1 2 3 pi/4 pi/6 pi/6 1 1]; % ====================================== options = optimset('Display','iter','MaxFunEvals',600,...     'TolFun',0.05,'TolX',0.05); %------------------- optim_param = fminsearch('altitude',param_init,options); % ====================================== save('paramametres_optim','optim_param');

The fminsearch options are set with optimset. Here, the maximum number of calls of the function altitude is 600. The optimal parameters are saved in the MAT-File parameters_optim at the end of the optimization process. We get, after several hours of computation:

optim_param = [1.0863,2.1026,3.2356,0.8892,0.4637,0.6572,0.9829,0.4884]

and indeed, running again bht_traject_compute

>> bht_traject_compute('datafilename','bird','controlsfilename','bird_cont','controlsparameters',optim_param);

and then bht_simulation:

>> film = bht_simulation('datafilename','bird','axis',[-11 11 -5 13],'displaytime','on','drawfishaxis',{3,'-r'});

we observe that the bird goes higher (11%) than it did previously (the resultat would have been more impressive over a larger time interval). Click on the image to play the movie: 

What about efficiency? (second optimization process)

Let us go back to our first flapping strategy. We wish now to improve the flapping efficiency rather than getting the bird to go upper. We begin with running the function bht_energy. Enter

>> bht_energy( 'datafilename','bird','timerange',[pi/2 6*pi-pi/2]);

to obtain the figure below. We set the option 'timerange' to [pi/2 6*pi-pi/2] to focus on the cruising regime.

A significantly part of the power expended by the bird is used to set the fluid (the surrounding air here) into motion. 

The quantity we want to lower is  

which reads, with our notation sum(sum(abs(torques)))*dT. The value of dT being fixed, we merely work with sum(sum(abs(torques))). For instance, when cont_parameters = [1 2 3 pi/4 pi/6 pi/6 1 1], we obtain:

sum(sum(abs(torques))) = 8.707938985714337*1e5.

An obvious way to reduce the torques should consist in slowing down the motion of the wings and even to not flapping at all. To avoid this problem, we impose the bird to go at least as high as it did with the initial value of cont_parameters. When  cont_parameters = [1 2 3 pi/4 pi/6 pi/6 1 1], we obtain:

Q(2,end)-Q(2,1) = 3.898381561990239.  

We modify the function altitude in order to obtain the following function torques_intensity:

function cost = torques_intensity(cont_parameters) % ========================================================================= ampl = cont_parameters(4:6); alpha(1:2) = cont_parameters(7:8); % test overlapping ======================================================== recov = 1.2+3.4*cos(ampl(1))+3.4*cos(ampl(1)+ampl(2)*(1+alpha(1)))...     +2.6*cos(ampl(1)+ampl(2)*(1+alpha(1))+ampl(3)*(1+alpha(2))); % computations ============================================================ if recov > 0.1     [T,p,dp,q,dq,Q] = bht_traject_compute('datafilename','bird',...         'controlsfilename','bird_cont',...         'controlsparameters',cont_parameters);     % ===================================     if Q(2,end)-Q(2,1)>=3.898381561990239     [T,ET,Efl,Efi,Epot,Edef,Enet,torques] = bht_energy('datafilename','bird',...             'drawgraph','off');     cost = sum(sum(abs(torques)))/1e5     else     %------------------------------     cost = 8.707938985714337;     %------------------------------     end     % =================================== else     %------------------------------     cost = 8.707938985714337;     %------------------------------ end

and we run the script:

% initial guess param_init = [1 2 3 pi/4 pi/6 pi/6 1 1]; % ====================================== options = optimset('Display','iter','MaxFunEvals',200,...     'TolFun',0.001,'TolX',0.001); %------------------- optim_param_torques = fminsearch('torques_intensity',param_init,options); % ====================================== save('parameters_optim_torques','optim_param_torques');

The maximum number of calls for the function torques_intensity is set to 200. Indeed, each evalutation of the function required to run both functions bht_traject_compute and bht_energy what means that the cost, in terms of CPU time, is quite high. After a couple of hours of computation, we get the result:

optim_param_torques = [0.7803 2.3271 3.1539 0.8015 0.5382 0.5699 0.9065 0.9838]

To see what this new flapping strategy looks like, enter:

>> bht_traject_compute('datafilename','bird','controlsfilename','bird_cont','controlsparameters',optim_param_torques);

and then run bht_simulation:

>> film = bht_simulation('datafilename','bird','axis',[-11 11 -5 13],'displaytime','on','drawfishaxis',{3,'-r'});

Click on the image to play the movie you obtain:

The bird reaches the same height as it did previously. Let us evaluate the efficiency of the flapping over the cruising regime:

>> bht_energy( 'datafilename','bird','timerange',[pi/2 6*pi-pi/2]);

The energy of the fluid is lower and so are the torques and the power expended at the joints.

Indeed, we obtain this time:

sum(sum(abs(torques))) = 7.01742031239794*1e5

which means a saving of 19.41% with respect to our initial flapping strategy.

2008 - A. Munnier and B. Pincon (Insitut Elie Cartan and INRIA Lorraine, Projet CORIDA, Nancy, France).