Code covered by the BSD License  

Highlights from
Solution of Differential Equations with MATLAB & Simulink: Lorenz Attractor Case Study

image thumbnail

Solution of Differential Equations with MATLAB & Simulink: Lorenz Attractor Case Study

by

 

11 Aug 2010 (Updated )

Simulink design pattern for solving differential equations, visualize results in MATLAB graphics

lorenz_graphs(block)
function lorenz_graphs(block)
% The goal of this code is to make the S-function behave like a Simulink
% block. During simulation of a model, the Simulink engine invokes this function at different times during the simulation
% Initialization step, update and output steps and the subfunctions
% corresponding to them get invoked in that same order. The initialization function gets invoked only once at the beginning of the 
% simulation. Update and output functions get invoked at every time step (major & minor). Also, observe that the level-2
% S-Function has a lot of similarity with C S-Function that the Level-1
% S-Function.
setup(block)

    function setup(block)
        % We setup and register callback methods that are directly called during simulation.
        
        % Number of input ports and dimensionality
        % (1) indicates the output port number-change it to add multiple ports.
        % Also do bookkeeping of these ports by documenting your code
        block.NumInputPorts  = 1;
        block.InputPort(1).Dimensions        = [1 3];
        block.InputPort(1).DirectFeedthrough = false;
        
        % Number of output ports and dimensionality
        block.NumOutputPorts = 0;
        % (1) indicates the output port number-change it to add multiple ports.
        % Also do bookkeeping of these ports by documenting your code
        % block.OutputPort(1).Dimensions= 0;
        
        % Number of block parameters
        block.NumDialogPrms = 1;
        block.DialogPrmsTunable ={'SimOnlyTunable'};
        
        % Number of continuous states, in our case we are assuming the state to be
        % discrete, so set this variable to zero. Note that this choice will impact
        % the choice of sample times. Try changing block.NumContStates = 1 and block.SampleTimes = [0 0] and you will
        % see that the plot looks very discrete as we are depending on the solver
        % to take the time steps. So, slowing the simulation down was a requirement
        % that could be enforced by making this block discrete and giving it a
        % sample time. The solver cannot bypass this "slow" block in the model. You
        % end up getting a nice plot with deterioration in simulation performance.
        block.NumContStates = 0;
        
        % Register the sample times.
        %  [0 offset]            : Continuous sample time
        %  [positive_num offset] : Discrete sample time
        %
        %  [-1, 0]               : Inherited sample time
        %  [-2, 0]               : Variable sample time
        block.SampleTimes = [0.01 0];
        
        %
        global tracer;
        
        % Register the various block methods with names of the functions mathcing
        % appropriately with the steps
        
        % Remember that DWork variables can only be set inside a
        % PostPropagationSetup method. So, we need to register that
        block.RegBlockMethod('PostPropagationSetup',@DoPostPropSetup);
        block.RegBlockMethod('InitializeConditions',@InitConditions);
        block.RegBlockMethod('Update', @Update);
        % Commenting this out as no output is required
        %block.RegBlockMethod('Outputs' ,@Output)
        
        % Code for SimState: idea of splitting simulations into stages and doing
        % that by using the end x of one stages as the x0 of the following stage. Look up the
        % documentation for more information on this
        block.SimStateCompliance = 'DefaultSimState';
        
        function DoPostPropSetup(block)
            % Setup Dwork
            % The PostPropagationSetup method initializes the DWork vector that stores the state
            % DWork vectors are blocks of memory that an S-function asks the Simulink engine to allocate to each instance of the S-function in a model.
            % If multiple instances of your S-function can occur in a model, your S-function must use DWork vectors instead of global or static memory to store
            % instance-specific values of S-function variables. Otherwise, your S-function runs the risk of one instance overwriting data needed by another instance,
            % causing a simulation to fail or produce incorrect results. The ability to keep track of multiple instances of an S-function is called reentrancy
            block.NumDworks = 1;
            % The (1) indicates this as the first state. You could add other states buy
            % subscripting (2), (3) and so on. Note that you have to do book-keeping
            % of these states so that you do not mix up-hence good to write some
            % documentation.
            block.Dwork(1).Name = 'x0';
            % DWork is a vector and so its dimension must be a positive integer
            block.Dwork(1).Dimensions      = 3;
            block.Dwork(1).DatatypeID      = 0;
            block.Dwork(1).Complexity      = 'Real';
            block.Dwork(1).UsedAsDiscState = false;
            
            
            function InitConditions(block)
                % Initialize Dwork
                block.Dwork(1).Data = (block.DialogPrm(1).Data)';
                % Visualization in a MATLAB GUI
                % Initialization code. Note that the GUI is being initialized as part of
                % the model execution. If you close the GUI during the Simulink model
                % execution, you will see a new GUI but this GUI will vanish because the
                % Simulink engine has no way of knowing that you undid what you set up
                % during the initialization step. Here is an enhanecment to think of:
                % Prevent the user from closing the GUI as long as the simulation is
                % running. How would you incorporate that?
                % Also, observe that the GUI code is not very sophisticated. The objects
                % that make this GUI are not connected in a rigid fashion to one figure
                % window. If there were multiple plotting functions inside your Simulink
                % model, this could cause a lot of overlap. Enhancement to think of: How to
                % design more robust GUIs so that they are delinked from each other?
                q=figure(1);
                set(q, 'Name', 'Level-2 Lorenz Oscillator Visualization', ...
                    'Color', [1 1 1]);
                set(gca,'Units', 'normalized');
                title('Position (x,y,z)');
                hold on;
                grid on
                axis([-50 50 -100 100 0 200])
                axis vis3d % freezes the aspect ratio to prevent auto strech to fill which can introduce jarring animation
                set(gca, 'xtick',[-50:25:50], 'ytick',[-100:50:100], 'ztick',[0:50:200]);
                set(gca,'xColor', [1 0 1], 'ycolor', [0 0 0], 'zcolor',[0 0 1])

                
                function Update(block)
                    % Connect the current computed point with the past connected point
                    line([block.InputPort(1).Data(1) block.Dwork(1).Data(1) ], ...
                        [block.InputPort(1).Data(2) block.Dwork(1).Data(2)], [block.InputPort(1).Data(3) block.Dwork(1).Data(3)],'LineWidth',1.5,...
                        'MarkerEdgeColor','b',...
                        'MarkerFaceColor',[.49 1 .63],...
                        'MarkerSize',12, 'EraseMode', 'none');
                    
                    % Plot an asterisk. At each update step, a new handle
                    % to the plot3 object is created. Use a simple delete 
                    % to create an effect similar to the comet3
                    % effect. See >>help COMET3 for more details.
                    a=plot3(block.InputPort(1).Data(1), block.InputPort(1).Data(2), block.InputPort(1).Data(3), 'r.', ...
                          'EraseMode', 'none', 'MarkerSize', 22);
                    delete(a);
                    
                    % Rotate the camera. Comment this code out if you do not want the orbiting
                    % effect
                    camorbit(0.5,0.5, 'data', [1 1 1]);
                    
                    % Update the points plotting
                    drawnow;
                    
                    % You want to hold the current input as the state(memory) for the next
                    % computation because you need it for the line command
                    block.Dwork(1).Data = block.InputPort(1).Data;
                    
                    % No output required here
                    %function Output(block)
                    % block.OutputPort(1).Data = block.Dwork(1).Data;

Contact us