How to apply A bang-bang signal of amplitude 1 N and 1 s width as an input force to reproduce certain results in Matlab?

Question

I working on dynamic modeling and simulation of a mechanical system (overhead crane), after I obtained the equation of motion, in the form: $$ M(q)\ddot{q}+C(q,\dot{q})\dot{q}+G(q)=Q $$

All the matrices are know inertia, $ M(q)$, Coriolis-Centrifugal matrix $ C(q,\dot{q})$, and gravity $ G(q)$ as functions of the generalized coordinates $q$, and their derivatives $\dot{q}$.

I want to solve for $q$, using Matlab ODE (in m-file), I got the response for some initial conditions and zero input, but, I want to find the response, for the aforementioned control signal (A bang-bang signal of amplitude 1 N and 1 s width), I'm trying to regenerate some results from the literature, and what the authors of that work said, regrading the input signal is the following: "A bang-bang signal of amplitude 1 N and 1 s width is used as an input force, applied at the cart of the gantry crane. A bang-bang force has a positive (acceleration) and negative (deceleration) period allowing the cart to, initially, accelerate and then decelerate and eventually stop at a target location." I didn't grasp what do they mean by bang-bang signal, I know in Matlab we could have step input, impulse, ...etc. But bang-bang signal, I'm not familiar with. According to this site and this bang bang is a controller rather.

Could anyone suggest to me how to figure out this issue and implement this input signal? preferably in m-file.

The code I'm using is given bellows, two parts:

function xdot = AlFagera(t,x,spec)
% xdot = zeros(8,1);
xdot = zeros(12,1); % to include the input torque

% % Crane Specifications
mp = spec(1);
mc = spec(2);
mr = spec(3);
L = spec(4);
J = spec(5);

g = 9.80;               % accelatrion of gravity (m/s^)

% % matix equations 
M11 = mr+mc+mp; M12 = 0; M13 = mp*L*cos(x(3))*sin(x(4)); M14 = mp*L*sin(x(3))*cos(x(4));
M21 = 0; M22 = mp+mc; M23 = mp*L*cos(x(3))*cos(x(4)); M24 = -mp*L*sin(x(3))*sin(x(4));
M31  = M13; M32 = M23; M33 = mp*L^2+J; M34 = 0;
M41 = M14; M42 = M24; M43 = 0; M44 = mp*L^2*(sin(x(3)))^2+J;
M = [M11 M12 M13 M14; M21 M22 M23 M24; M31 M32 M33 M34; M41 M42 M43 M44];

C11 = 0; C12 = 0; C13 = -mp*L*sin(x(3))*sin(x(4))*x(7)+mp*L*cos(x(3))*cos(x(4))*x(8);
C14 = mp*L*cos(x(3))*cos(x(4))*x(7)-mp*L*sin(x(3))*sin(x(4))*x(8);
C21 = 0; C22 = 0; C23 = -mp*L*sin(x(3))*cos(x(4))*x(7)-mp*L*cos(x(3))*sin(x(4))*x(8);
C24 = -mp*L*cos(x(3))*sin(x(4))*x(7)-mp*L*sin(x(3))*cos(x(4))*x(8); 
C31 = 0; C32 = 0; C33 = 0; C34 = -mp*L^2*sin(x(3))*cos(x(3))*x(8);
C41 = 0; C42 = 0; C43 = -C34; C44 = mp*L^2*sin(x(3))*cos(x(4))*x(7);

C = [C11 C12 C13 C14; C21 C22 C23 C24; C31 C32 C33 C34; C41 C42 C43 C44];
Cf = C*[x(5); x(6); x(7); x(8)];

G = [0; 0; mp*g*L*sin(x(3)); 0];

fx = 0; 

if t >=1 && t<=2
fy = 1.*square(t*pi*2);
else fy = 0;
end

F =[fx; fy; 0; 0];     % input torque vector, 

xdot(1:4,1)= x(5:8);
xdot(5:8,1)= M\(F-G-Cf);
xdot(9:12,1) = F;

And:

clear all; close all; clc;

t0 = 0;tf = 20;

x0 = [0.12 0.5 0 0, 0 0 0 0,0 0 0 0];  % initional conditions

% % spectifications
Mp = [0.1 0.5 1];      % variable mass for the payload
figure
plotStyle = {'b-','k','r'};
for i = 1:3
mp = Mp(i);
mc = 1.06; mr = 6.4;           % each mass in kg
L = 0.7; J = 0.005;            % m, kg-m^2 respe.
spec = [mp mc mr L J];
% % Call the the function
[t,x] = ode45(@(t,x)AlFagera(t,x,spec),[t0 :0.001: tf],x0);

legendInfo{i} = ['mp=',num2str(Mp(i)),'kg'];


fx = diff(x(:,9))./diff(t);
fy = diff(x(:,10))./diff(t);
tt=0:(t(end)/(length(fx)-1)):t(end); % this time vector

% to plot the cart positions in x and y direcitons
subplot(1,2,1)
plot(t,x(:,1),plotStyle{i})
axis([0 20 0 0.18]);
grid
xlabel('time (s)');
ylabel('cart position in x direction (m)');
hold on
legend(legendInfo,'Location','northeast')

subplot(1,2,2)
plot(t,x(:,2),plotStyle{i})
axis([0 20 0 1.1]);
grid
xlabel('time (s)');
ylabel('cart position in y direction (m)');
hold on
legend(legendInfo,'Location','northeast')

end

% to plot the input torque, (bagn-bang signal), just one sample
figure
plot(tt,fy)
grid
set(gca,'XTick',[0:20])
xlabel('time (s)');
ylabel('input signal, f_y (N)');

Furthermore, the results I'm getting and what I supposed to get are shown:

Major difficulties, initial conditions are not clearly stated in the paper, the input force direction, is only in y (which it should be), or it has different direction. I appreciate any help.

the paper I'm trying to recreate is: R. M. T. Raja Ismail, M. A. Ahmad, M. S. Ramli, and F. R. M. Rashidi, “Nonlinear Dynamic Modelling and Analysis of a 3-D Overhead Gantry Crane System with System Parameters Variation.,” International Journal of Simulation–Systems, Science & Technology, vol. 11, no. 2, 2010. http://ijssst.info/Vol-11/No-2/paper2.pdf

Answer 1

This is how I would go about simulating a nonlinear ODE in Matlab. As I mentioned in a (now-deleted) comment on your question, I typically work with linear ODE's in my line of work, which means that I usually use the (awesome) functions in the Control System Toolbox.

Now, I'll start by saying that you haven't given any definitions of what the matrices are in your equation, so I can't give you a concrete example of how to solve your specific example. If you update your question to give numbers for everything then I can modify this answer to solve your problem. That said, you give:

$$ M(q)\ddot{q}+C(q,\dot{q})\dot{q}+G(q)=Q $$

First I would say that this can be re-written as:

$$ \ddot{q} = \frac{1}{M(q)}\left(Q-C(q,\dot{q})\dot{q} -G(q)\right) \\ $$

So, for an example for you, I chose an RLC circuit, which takes the form:

$$ \begin{array} .\dot{I} = \frac{1}{L}\left(V_{\mbox{in}} - V_C - IR\right) \\ \dot{V} = \frac{1}{C}I \end{array} $$

Typically your input signal would be a smooth function of time. Here you're looking for a bang-bang signal, which is akin to a light switch. At some time, the input signal goes from nothing immediately to some value, then later from that value immediately back to nothing.

So, where typically you would use an interpolate command to get values defined sample time increments, here you don't want those commands interpolated. You don't want a ramp, you want the signal to immediately shift.

So this question really is two parts:

1. How do I pass parameters to an ODE45 (or 23 or whatever else) function in Matlab, and
2. How do I define a step change in an input signal for an ODE45 function?

The answers are (examples to follow)

1. Define your ODE function as a function in its own script, such that the first line of the script is something like function [dq] = AlFageraODE(t,x,parameter1,parameter2,...,parameterN). Then, when you want to solve the ODE, call Matlab's built-in ODE solver as follows: [t,y] = ode45(@AlFageraODE,[t0,t1],[y0,y1],[ ],parameter1,parameter2,...,parameterN);. The square brackets in bold there need to be included because that is where you can pass specific options to the ODE solver. If you don't want to pass anything specific and are okay with the default settings, you still need to put something, so put an empty array - this is done with an empty array []. After that you can put in parameters that will be passed to your custom function.
2. To get a true step function, you need to split the simulation into three distinct sets - before the step, during the step, and after the step. Anything else will result in the need to interpolate the input command. The last entry in the outputs will be the initial conditions for the next segment of the simulation.

Below, I've written the custom ODE function, and then below that is the script used to solve the custom ODE. Note that the R, L, and C values are set in the calling script, not in the custom ODE. This is because those values are passed in to the custom ODE, along with what the applied voltage should be during that particular segment.

The examples:

1. The custom ODE function (again, this simulates an RLC circuit, but you can modify it for any custom ODE needing parameters passed in.)

.

function dx = RoboticsODE(t,x,Vin,R,L,C)
dx = zeros(2,1);

I = x(1);
V = x(2);

di_dt = (1/L)*(Vin - V - I*R);
dv_dt = (1/C)*I;

dx(1) = di_dt;
dx(2) = dv_dt;

2. The script that solves that ODE for a bang-bang signal.

.

%% Clear/close everything
clear;
close all;
clc;

%% Simulation Configuration
v0 = 0;             % Initial *output* voltage
i0 = 0;             % Initial circuit current
C = 47*(1/1000);    % Capacitance, F
R = 0.1;            % Resistance, Ohm
L = 22*(1/1000);    % Inductance, H

appliedVoltages = [0,1,0]; % The applied voltages at t0, bangOn, bangOff
tStart = 0;         % t0
bangStart = 1;      % bangOn time
bangWidth = 1;      % bangOff time = bangOn time + bangWindow
endWindow = 5;      % endWindow is how long to "watch" the simulation after 
                    %   the bang-bang signal goes "off".

%% Output Initialization
outputTime = zeros(0,1);
outputVoltage = zeros(0,1);
outputCurrent = zeros(0,1);
inputVoltage = zeros(0,1);

%% Dependent Configuration
currentValues = [i0;v0];
samplePoints = cumsum([tStart,bangStart,bangWidth,endWindow]);
% A note on the above - I use the cumulative sum (cumsum) because I defined
% each point as the number of seconds *after* the previous event. If you
% defined absolute time points then you'd just use those time points
% directly with no cumulative sum.
nSegments = numel(samplePoints)-1;

%% Simulation
for currentSegment = 1:nSegments
    % Setup the simulation by getting the current time window and "intial"
    % conditions.
    Vt = appliedVoltages(currentSegment);
    t0 = samplePoints(currentSegment);
    t1 = samplePoints(currentSegment+1);
    sampleTime = [t0;t1];

    % Run the simulation by solving the ODE for this particular segment.
    [intermediateTime,intermediateOutput] = ...
        ode45(@RoboticsODE,sampleTime,currentValues,[],Vt,R,L,C);

    % Assign outputs
    nOutputPoints = numel(intermediateTime);
    outputTime(end+1:end+nOutputPoints) = intermediateTime;
    outputCurrent(end+1:end+nOutputPoints) = intermediateOutput(:,1);
    outputVoltage(end+1:end+nOutputPoints) = intermediateOutput(:,2);
    inputVoltage(end+1:end+nOutputPoints) = Vt*ones(nOutputPoints,1);

    % Setup the next simulation by setting the "initial" conditions for
    % that simulation equal to the ending conditions for the current
    % simulation.
    currentValues(1) = outputCurrent(end);
    currentValues(2) = outputVoltage(end);
end

%% Output Plot
plot(outputTime,inputVoltage);
hold on;
plot(outputTime,outputVoltage);
title('RLC Circuit with Step Input to ODE45');
xlabel('Time (s)');
ylabel('Voltage (V)');
legend('Input Voltage','Output Voltage');

3. The plot of the output.

Finally, as I mentioned, if you would be willing to give concrete numbers for your equation then I could tailor this to solve your particular problem. As it stands, I can't provide any solution to a symbolic nonlinear ODE so this example is the best I can give you.

:EDIT:

I've got the problem solved for you. The code is attached below. I'll say this, though: It's important, for a step input (bang-bang, etc.) to segment the ODE solver process, as I described above. This is because Matlab tries to optimize solving the ODE and may not take process time in exactly the way you would expect.

The segmenting method is again, as described above, where you split the solving process at every discontinuity. The initial conditions for the following step is equal to the ending conditions of the current step.

The images below are the solutions I got with the segmented method and the "all-in-one" method. The all-in-one is the way you had it setup, where the force function was determined by an if statement in the ODE function. Because Matlab chooses the sample time increments, the positive and negative segments aren't guaranteed to have exactly the same number of samples. I think this is the reason for the drift in the output of the all-in-one solution.

I found several problems with your method. When I corrected them, I got a plot that looked (to me) to exactly duplicate the plots from the paper you linked.

1. The biggest problem - fx and fy should be the same.
2. Also a problem, the pulse width should be 1s. You used a square wave of period 1s, meaning that there was a 0.5s "positive" and 0.5s "negative" signal. I halved the frequency and got the proper 1s signal width.
3. Your initial conditions for x and y were not zero. They are zero in the paper, so I set them to zero in the simulation to replicate the figures in the paper.

Here's the code! First, the ODE script:

.

function varDot = AlFagera(t,var,spec,F)
% In general, I did a lot of cleanup with this function to make things
% easier for me to read.

%% Misc Declarations
% varDot = zeros(8,1);
varDot = zeros(12,1); % to include the input torque
g = 9.80;               % Acceleration of gravity (m/s^2)

%% Define Forces if Undefined
% If the segmentSolution is being used, then the force is supplied to the
% function. If it's not being used (the "all-in-one" solution), then the
% force is determined by the simulation time. 
if isempty(F)
    if t >=1 && t<=3
        fy = -1.*sign(sin(t*pi));
    else
        fy = 0;
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % I'll highlight this, because I think this was the key problem     %
    % you were having. The force for fx and fy should be the same. You  %
    % had fx = 0 for all cases.                                         %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %fx = 0;
    fx = fy;
    F = [fx;fy;0;0];
end
%% Crane Specifications
mp = spec(1);
mc = spec(2);
mr = spec(3);
L = spec(4);
J = spec(5);

%% Breakout the Input Variable
x        = var(1);
y        = var(2);
theta    = var(3);
phi      = var(4);

xDot     = var(5);
yDot     = var(6);
thetaDot = var(7);
phiDot   = var(8);

%% Repackage the Inputs into Useful Form
q = [...
    x; ...
    y; ...
    theta; ...
    phi];

qDot = [...
    xDot; ...
    yDot; ...
    thetaDot; ...
    phiDot];

%% Simplified expressions for (extensive) use later
mpL  = mp*L;
cosT = cos(theta);
sinT = sin(theta);
cosP = cos(phi);
sinP = sin(phi);

%% Matrix Expressions
M = [...
    mr+mc+mp      , 0              , mpL*cosT*sinP , mpL*sinT*cosP; ...
    0             , mp+mc          , mpL*cosT*cosP , -mpL*sinT*sinP; ...
    mpL*cosT*sinP , mpL*cosT*cosP  , mpL*L+J       , 0; ...
    mpL*sinT*cosP , -mpL*sinT*sinP , 0             , mp*(L*sinT)^2+J];
C = [...
    0 , 0 , -mpL*sinT*sinP*thetaDot+mpL*cosT*cosP*phiDot , mpL*cosT*cosP*thetaDot-mpL*sinT*sinP*phiDot; ...
    0 , 0 , -mpL*sinT*cosP*thetaDot-mpL*cosT*sinP*phiDot , -mpL*cosT*sinP*thetaDot-mpL*sinT*cosP*phiDot; ...
    0 , 0 , 0                                            , -mpL*L*sinT*cosT*phiDot; ...
    0 , 0 , mpL*L*sinT*cosT*phiDot                       , mpL*L*sinT*cosT*thetaDot];
G = [...
    0; ...
    0; ...
    mpL*g*sinT; ....
    0];

%% Assign Outputs
qDdot = M\(-C*qDot-G+F);
varDot = [...
    qDot; ...
    qDdot; ...
    F];

Then, the script that solves the ODE:

.

clear all; 
close all; 
clc;

% Compare the all-in-one method of solving the problem with the segmented
% method of solving the problem by setting the variable below equal to
% "true" or "false". 
segmentSolution = true;

t0 = 0;tf = 20;

% Your initial conditions here for x- and y-positions were not zero, so I
% set them to zero to reproduce Figure 2 and Figure 3 in the paper you
% linked. 
% Also - you don't ever use the last 4 values of this except as a way to
% output the force. This isn't used in the segmentSolution because there
% the input force is supplied to the function. 

x0 = [0 0 0 0, 0 0 0 0,0 0 0 0];  % Initial Conditions

%% Specifications
Mp = [0.1 0.5 1];      % Variable mass for the payload
figure
plotStyle = {'b-','k','r'};

%% SegmentSolution Settings
fx = [0,1,-1,0];
fy = [0,1,-1,0];
tStart = 0;
tOn = 1;
bangWidth = 1;
tEndWindow = 17;
sampleTime = cumsum([tStart,tOn,bangWidth,bangWidth,tEndWindow]);
nSegments = numel(sampleTime)-1;
%% Simulation
for i = 1:3
    mp = Mp(i);
    mc = 1.06; mr = 6.4;           % each mass in kg
    L = 0.7; J = 0.005;            % m, kg-m^2 respe.
    spec = [mp mc mr L J];
    %% Call the the function
    initialConditions = x0;

    if segmentSolution
        t = zeros(0,1);
        x = zeros(0,numel(x0));
        outputFx = zeros(0,1);
        outputFy = zeros(0,1);
        for currentSegment = 1:nSegments
            inputForce = [fx(currentSegment),fy(currentSegment),0,0].';
            t0 = sampleTime(currentSegment);
            t1 = sampleTime(currentSegment+1);
            [intermediateT,intermediateX] = ode45(@AlFagera,[t0 :0.001: t1],initialConditions,[],spec,inputForce);
            nOutputSamples = numel(intermediateT);
            index1 = size(t,1)+1;
            index2 = size(t,1)+nOutputSamples;
            t(index1:index2) = intermediateT;
            x(index1:index2,:) = intermediateX;
            initialConditions = x(end,:).';
            outputFx(index1:index2) = inputForce(1)*ones(nOutputSamples,1);
            outputFy(index1:index2) = inputForce(2)*ones(nOutputSamples,1);
        end
        tt = t;
    else
        inputForce = []; % Leave this empty for the all-in-one solver. 
                         % There's a check in the code to setup the force
                         % when it's not defined. 
        [t,x] = ode45(@AlFagera,[t0:0.001:tf],initialConditions,[],spec,inputForce);
        outputFx = diff(x(:,9))./diff(t);
        outputFy = diff(x(:,10))./diff(t);
        tt=0:(t(end)/(length(outputFx)-1)):t(end);
    end
    legendInfo{i} = ['mp=',num2str(Mp(i)),'kg'];


    %fx = diff(x(:,9))./diff(t);
    %fy = diff(x(:,10))./diff(t);
    %tt=0:(t(end)/(length(fx)-1)):t(end); % this time vector

    % to plot the cart positions in x and y direcitons
    subplot(1,2,1)
    plot(t,x(:,1),plotStyle{i})
    axis([0 20 0 0.18]);
    grid
    xlabel('time (s)');
    ylabel('cart position in x direction (m)');
    hold on
    legend(legendInfo,'Location','northeast')

    subplot(1,2,2)
    plot(t,x(:,2),plotStyle{i})
    axis([0 20 0 1.1]);
    grid
    xlabel('time (s)');
    ylabel('cart position in y direction (m)');
    hold on
    legend(legendInfo,'Location','northeast')
end

% to plot the input torque, (bang-bang signal), just one sample
figure
plot(tt,outputFy);
grid
set(gca,'XTick',[0:20])
xlabel('time (s)');
ylabel('input signal, f_y (N)');

I took the liberty of cleaning up the code in the ODE function, I hope you don't mind.

So, to summarize, the solution to your problem is:

1. Pass inputs to the ODE function (such as parameters and applied forces), and
2. Segment the ODE solution at each discontinuity.

Answer 2

In bang-bang control, control inputs can only be activated similar to a switch: when they are engaged they are fully on positive, or fully on negative. Some bang-bang control studies (not many that I can recall) account for a slew rate of the controller's output, so the outputs are not modelled as square waves, but rather as trapezoids. I think this is the "bang-bang signal" you need to generate. The phrasing is odd.

It sounds as if you need to add a control force to your model that is just a trapezoid (acceleration period, constant full-on period, and deceleration period). In m-code it would be something like the following. Here I am assuming your step size is 1 msec, with a 10% slew time. You may want different slew parameters. I've written it so you can change step size and slew percents in one place.

max_F = 1.0;
steps_per_sec = 1000;
slew_up = 0.1;
slew_down = 0.1;
accel_steps = floor(steps_per_sec * slew_up);
decel_steps = floor(steps_per_sec * slew_down);
full_on_steps = steps_per_sec - accel_steps - decel_steps;

loops = 1;
// acceleration period
for i = 1:accel_steps
    F(i) = i * max_F / accel_steps;
    loops = loops + 1;
end
// full-on period
for i = loops : accel_steps + full_on_steps
    F(loops) = max_F;
    loops = loops + 1;
end
// deceleration period
for i = loops : accel_steps + full_on_steps + decel_steps
    F(loops) = max_F - (i - accel_steps - full_on_steps) * max_F / decel_steps;
    loops = loops + 1;
end
// after input is back to zero
for i = loops : 10 * steps_per_sec
    F(i) = 0.0;
end

// now apply F to your simulation

@AlFagera, Thanks for posting a reference to the actual paper. Following the references and re-searching led me to this paper: M. S. Alam, M. A. Hossain, and M. O. Tokhi, Design of Command Shaper Using Gain-Delay Units and Particle Swarm Optimisation for Vibration Control of a Flexible System, 14th International Congress on Sound & Vibration, Cairns, Australia 2007 (https://www.acoustics.asn.au/conference_proceedings/ICSV14/papers/p197.pdf). In that paper they describe a bang-bang signal as I've done above, except their acceleration and deceleration periods are quite short.

Answer 3

The 'bang' refers to a step transition. Bang-Bang means that one applies maximum positive effort, then maximum negative effort, with step transitions in between.

Technically we starting at 0 velocity and returning to 0 velocity so it should be bang-bang-bang control; but it sounds strange to say that.

For theoretical analysis involving bang-bang inputs, I would use the dirac delta. I don't use matlab that much but it looks like they have it for symbolic analysis and numerical analysis. It sounds like you already know this though.