Synthesizing a feedforward control stragegy for stiff hybrid systems

This demo shows how a model of a variable-step solver can be exploited to make stiff hybrid dynamic systems amenable to model checking technology. Specifically, a feedforward control strategy will be synthesized for an operation by a surface mount device that is numerically stiff because of widely different time scales.

Details about the theoretical background can be found in the article "A Computational Model of Time for Stiff Hybrid Systems Applied to Control Synthesis" by P.J. Mosterman, J. Zander, G. Hamon, and B. Denckla as published in the August 2011 issue of the IFAC journal Control Engineering Practice.


The surface mount device system

The overall system comprises subsystems for the numerical variable-step solver (cyan), the control and surface mount device system dynamics (blue), the control properties (purple), and the data visualization (gray).

smd_top = 'surface_mount_device_R11a_v1';
smd_dv = 'surface_mount_device_dv_R11a_v1';
smd_harness = 'surface_mount_device_dv_R11a_v1_harness';


A model of the surface mount device

A surface mount device mounts a component on a board in three stages. It first quickly moves to make contact with the board. In the second stage, it holds the component against the board with low acceleration. This holding must last for a minimum duration after which in the third stage the component is released.


The physics are modeled as a component with a given mass upon which two forces are acting: (i) a control force and (ii) a contact force.

open_system([smd_top, '/surface mount device']);

The contact force is only active when the component that is mounted is in contact with the board. When in contact the force is derived from a stiff spring/damper system. The contact behavior introduces a discontinuity when switching from the mode where there is no contact to the mode where there is contact. To properly effect this mode switch, a rate transition synchronizes the contact determination with the numerical solver by means of the 'contact' and 'contact_rt' port pair.

open_system([smd_top, '/surface mount device/contact behavior']);

hilite_system([smd_top, ...
    '/surface mount device/contact behavior/contact']);
hilite_system([smd_top, ...
    '/surface mount device/contact behavior/contact_rt']);

The control force

The force that controls the stamper motion of the surface mount device is generated as a feedforward force profile. This profile is time based and because of applying the explicit model of a variable-step solver, time is computed as an explicit variable in the model. This explicit computation allows a latch block to maintain the time of the most recent accepted integration step as well as the corresponding position.

open_system([smd_top, '/control']);

In more detail, the control force also requires the time at which contact with the board is made, as that is when a precise sequence of forces is exerted for specific respective durations.

open_system([smd_top, '/control/control force profile']);

The sequence of forces is produced by moving through a series of controller states as captured by a Stateflow chart. The transition from one control force to the next is made based on the necessary duration for each of the control forces. Because of the use of absolute time, the enabling conditions contain the cumulative duration at each state in the sequence.

open_system([smd_top, '/control/control force profile/control force']);

Computing the behavior

To have control over the step in time that is made, the solver computes time explicitly by integrating a constant input of value 1. Because time is computed by the variable-step solver subsystem, the Simulink solver is set to FixedStepDiscrete with a nominal value of 1.

Note that none of the blocks that are used in the model employ time as computed by the internal Simulink solver. Indeed, the independent variable of Simulink should now not be interpreted as time but, rather, as the number of evaluations in the overal computation of the behavior. Time is represented differently and the details can be found in the suggested reading.

This is why, a block such as the Ramp source block cannot be used because it computes its output value based on the current time as determined by the internal Simulink solver and not the explicitly modeled variable-step solver.

The variable-step solver

Because of the spring/damper system that becomes active when the component to be mounted makes contact with the board, the dynamics of the overall behavior change from relatively slow while approaching the board to relatively fast when in contact with the board. Systems with such widely varying time constants in their dynamics may be considered 'stiff' from a numerical perspective. In other words, numerical integration must proceed with relatively small steps during the fast behavior. A variable-step solver dynamically adapts the step size of the numerical integration to solve such stiff systems efficiently.

The variable-step solver subsystem contains rate transition behavior that is necessary to properly combine continuous-time behavior with discrete-time behavior.

open_system([smd_top, '/variable step solver']);

Variable-step integration

In general, there are many approaches to variable-step integration. The variable-step integration employed here computes the integrated values over an integration time step in two stages. First, a Forward Euler approximation of the values is computed by the compute Euler approximation subsystem.

open_system([smd_top, '/variable step solver/variable-step integration']);

The Forward Euler approximation computes the new values of the state x over a time step of size h by moving along the gradient (shown in blue) at the current time. This results in an unknown error.

title('Forward Euler Approximation','fontsize',14);

Second, a more accurate Trapezoidal-based approximation is computed by the compute Trapezoidal approximation subsystem. The Trapezoidal-based approximation takes the average of the gradient at the current time (shown in blue) and after one time step (shown in green). Because the values after one time step are unknown, the gradient in green is computed based on the Forward Euler approximation and this (explicit) implementation of a Trapezoidal-based approximation is also referred to as a Heun approximation.

title('Trapezoidal Approximation','fontsize',14);

The new values of x are then computed by following the average of the two gradients (shown in red) over one time step, h. This results in another error that is unknown.

The Forward Euler and Trapezoidal increments are then compared by the compute step size subsystem. If there is a discrepancy between them that exceeds a given tolerance, the step size is reduced and the computations are repeated to compute an alternate and more accurate state.

The output of the solver alternates between the Euler approximation and the Trapezoidal approximation, as implemented by the alternate between Euler and Trapezoidal subsystem. This allows the Trapezoidal approximation to compute the necessary end point gradient based on the Euler approximation.

When the discrepancy falls below the given tolerance, the integration point is accepted and the corresponding time is said to be a major time step (as opposed to a minor time step). Time can now be advanced and the next integration point is computed.

Executing the system

Now, the overall system behavior can be analyzed by starting the simulation. First, the position of the component at each of the evalations by the variable-step solver can be studied. For clarity, the data points in red circles are the state values as computed within the prescribed integration tolerance of the variable-step solver (the accepted integration values). It can be seen that around the point of contact (at 0), the required tolerance necessitates many evaluations where the computed step in the position value is repeatedly reduced till the accepted values are close to 0.


plot(td(:,1), xd(:,2), 'b*');
grid on;
hold on;
plot(tm(:,1), xm(:,2), 'ro');
title('Position vs. Evaluations','fontsize',14)
hold off;

The advance of time against the number of evaluations shows that when the component is not subject to the large board force, time advances quickly with a constant rate. When the contact behavior is active, time advances relatively slowly. Zooming in on when the contact behavior becomes active, it can be seen that time recedes repeatedly in an attempt to achieve the desired accuracy as specified by the tolerance of the solver.

Also notice that time is kept constant for pairs of consecutive evaluations. This reflects the two-stage nature of the solver, which implements first the Euler approximation and then the Trapezoidal approximation, before the discrepancy between the two can be evaluated.

plot(td(:,1), td(:,2), 'b*');
grid on;
title('Progress of Time','fontsize',14)

The ultimate objective of the control force profile is to keep the acceleration low for a specific duration when the component makes contact with the board.

plot(td(:,2), ad(:,2), 'b*');
grid on;
hold on;
plot(tm(:,2), am(:,2), 'ro');
title('Acceleration vs. Time','fontsize',14)
hold off;

Of particular interest here are the times at which the computed state values are accepted by the variable-step solver as satisfying the tolerance requirement.

When zooming in on the contact behavior, it can be seen that the first time step for which the absolute acceleration is kept below 0.75 starts at tm(11,2) which is time 0.031875. The absolute acceleration next exceeds 0.75 at tm(16,2) which is time 0.038125. The difference of 0.00625 satisfies the required minimum duration of 6.25 ms.

grid on; hold on;
title('Accepted Acceleration vs. Time in Detail','fontsize',14)
axis([0.03 0.04 -20 5]);

hMaxAccel = plot([0.03 0.045],[0.75 0.75],'g-.','LineWidth',2);
hMinAccel = plot([0.03 0.045],[-0.75 -0.75],'g-.','LineWidth',2);
hAccelGroup = hggroup;

hMaxDuration = plot([tm(16,2) tm(16,2)], [-20 5],'b-.','LineWidth',2);
hMinDuration = plot([tm(11,2) tm(11,2)], [-20 5],'b-.','LineWidth',2);
hStairsDuration = plot(...
    [tm(11,2) tm(11,2) tm(12,2) tm(12,2) tm(13,2) tm(13,2) tm(14,2) ...
        tm(14,2) tm(15,2) tm(15,2) tm(16,2)], ...
    [am(11,2) am(12,2) am(12,2) am(13,2) am(13,2) am(14,2) am(14,2) ...
        am(15,2) am(15,2) am(16,2) am(16,2)], ...
hDurationGroup = hggroup;

plot(tm(:,2), am(:,2), 'ro','MarkerFaceColor','w','LineWidth',2);

legend('acceleration limits',...
hold off;

Synthesizing the control profile

The control profile that generated this particular acceleration shows a force of -10 during the first phase. In the second phase, a sequence of specific control forces with a value of either -3.5, -4, or -4.5 is exerted. Finally, in the third phase, when the component is released, the stamper is withdrawn with a control force of 10.

plot(tm(:,2), fm(:,2), 'ro','MarkerFaceColor','w','LineWidth',2);
grid on;
hold on;
title('Control Force vs. Time','fontsize',14)
ylabel('control force');
axis([0.03 0.04 -10 10]);
hold off;

Simulink Design Verifier model

The sequence of control forces in the second phase must be carefully determined. Here, a 'model checking' approach is employed by using Simulink Design Verifier to prove that the component cannot be held against the board with a specific low acceleration for a minimum duration. Simulink Design Verifier then attempts to prove this property. If it can falsify the property, Simulink Design Verifier produces a counterexample. The force profile of the counterexample is exactly the sought after force sequence.

To create a model for Simulink Design Verifier, the control force that can be explored to falsify the property is assumed to be one of either -3.5, -4, and -4.5 at any time.


The complexity of the model checking problem is reduced by considering the point at which the component starts being pressed against the board. This allows the mode switch to the contact model to be set to always true.

The initialization of the model must then employ the values computed for the point of contact. Inspect xm of the previous simulation and notice how the first value below 0 is at index 11, which is evaluation 48. So, to initialize use the time, position, and velocity values at index 11.

format long e
xm =

  Column 1


  Column 2


ans =


ans =


ans =


Running Simulink Design Verifier

Now start proving the property by running Simulink Design Verifier from the Tools menu and select 'Prove Properties'. This will falsify the property and create a test harness to generate a counter example. Instead of performing the verification, which may take a long time (depending on your hardware, a couple of hours), you can open up the pregenerated test harness model.


Translating the control profile into a Stateflow chart

To bring the input that generates the desired control profile back into the original model, compute the durations of each of the control force values. Run the test harness. By looking at fm, you can tell when the force changes and then compute the duration by determining the difference from the previous time. The cumulative durations then become transition conditions for the Stateflow chart that produces the control force profile in the original model.

simOut = sim(smd_harness,'StopTime','36');

fm = simOut.get('fm');
tm = simOut.get('tm');
format short e;

tm(2,2) - tm(1,2)

tm(4,2) - tm(2,2)

tm(5,2) - tm(4,2)

tm(7,2) - tm(5,2)
ans =


ans =


ans =


ans =


ans =


ans =


ans =


ans =