Document Sample
					                                                                                                                               2008 ANNUAL REPORT           19

                                                              A. INDUSTRIAL CONTROL
                                                                 AND OPTIMISATION
                                                              Program Goals:
                                                              The partnerships between researchers and industry enable reciprocal transfer of
                                                              knowledge and new ideas of great potential impact on the community and economy.
                                                              This program encompasses several research projects motivated by and in collaboration
                                                              with industrial partners. The main underlying theme of these projects is the application
                                                              of advanced modelling, control and optimisation techniques to maximise asset utilisation
                                                              and improve performance. The complexity of the dynamics of such processes arise from
                                                              factors including model errors, unknown disturbances, nonlinearities, distributed parameter
                                                              systems, elements of Human-Machine Interaction and hybrid (discrete and continuous
                                                              state) components. Expected outcomes of the program include high quality research
                         Julio Braslavsky                     solutions and human resources tailored to the needs of Australian industry.
                         Program Leader
                                                              In 2008 we welcomed CfW Hamilton Jet & Co (New Zealand) and Boeing Research &
                                                              Technology, Australia as Industrial Affiliates, joining Halcyon International, CSR, Connell
                                                              Wagner and Industrial Automation Services (Hatch IAS).

                                                              A.1 PERfORMANCE OPTIMISATION Of MARINE SYSTEMS
                                                              Project Leader: T. Perez
                                                              Researchers:      C. Løvaas, G.C. Goodwin and J.C. Agüero
                                                              External Academic Collaborators:
                                                              T.I. fossen (Norwegian University of Science and Technology, Norway)
                                                              C. Holden (Norwegian University of Science and Technology, Norway)
                                                              E. Herrero, (Student, University of Cantabria, Spain)
                         Tristan Perez
                         Deputy Program Leader                External Industrial Collaborators:
                                                              P. Steinmann (Halcyon International, Australia)
                                                              T. Armstrong (Austal Ships, Australia)
                                                              J. Borret (Hamilton Jet, New Zealand)
                                                              M. Santos-Mujica (Robotiker-Tecnalia, Spain)
                                                              S. Peder-Berge (Offshore Simulator Centre, Norway)
                                                              Marine systems are designed to perform complex operations that require appropriate
                                                              reliability and economy. These requirements demand an interdisciplinary approach to
                                                              address the tight integration of design aspects related to hydrodynamics, structures, and
                                                              motion control.
                                                              This project is dedicated to the design of tools for guidance and motion control with the
                                                              aim of optimising the performance of marine vehicles in different operations. The project
                                                              targets vessels and operations within offshore, maritime transport, underwater exploration,
                                                              unmanned vehicles and wave energy conversion. Some of the current research is being
                                                              conducted together with international academic and industry collaborators.

 figure 1: Marine systems are designed to perform
        complex operations that require appropriate
        reliability and economy. These requirements
   demand an interdisciplinary approach to address
   the tight integration of design aspects related to
     hydrodynamics, structures, and motion control.
(Picture courtesy of Offshore Simulator Centre AS, Norway).
20   ARC Centre of Excellence for Complex Dynamic Systems and Control

A                                                    A.1.1 Marine system simulation tools
                                                     Researchers: T. Perez and T. fossen (Norway)
                                                     The Marine Systems Simulator (MSS) is an environment developed to provide the
                                                     necessary resources for rapid implementation of mathematical models of marine systems
                                                     with focus on control system design. The platform adopted for the development of the
                                                     MSS is Matlab/ Simulink. This allows a modular simulator structure, and the possibility of
                                                     distributed development. Openness and modularity of software components have been
                                                     prioritised in the design. This enables a systematic reuse of knowledge and results in
                                                     efficient tools for research and education.
                                                     In 2008, we added functions for frequency-domain identification, which are an aid to
                                                     the construction of time-domain models from frequency-domain data computed from
                                                     hydrodynamic codes. The latest version of the software and future updates can be freely
                                                     downloaded from

                                                     A.1.2 Identification of radiation force parametric models of
                                                     marine structures from 2D frequency domain data
                                                     Researchers: T. Perez and T. fossen (Norway)
                                                     The ability to predict ship responses and loads in waves is an important tool in the
                                                     design of marine structures and motion control systems. One method for constructing
                                                     time-domain models consists of using the data generated by the hydrodynamic codes
                                                     to compute the different elements of the so called Cummins’ equation of ship motion. In
                                                     this project, we have been studying the application of different identification methods in
                                                     both time and frequency domain to make best use of the available hydrodynamic data and
                                                     constraints based on prior knowledge derived from hydrodynamic theory.
                                                     In 2008, we addressed the problem of joint identification of infinite-frequency added mass
                                                     and fluid-memory models of marine structures from finite frequency data. This problem
                                                     is relevant for cases where the code used to compute the hydrodynamic coefficients of
                                                     the marine structure does not give the infinite frequency added mass. This case is typical
                                                     of codes based on 2D-potential theory codes. The method proposed presents a simpler
                                                     alternative approach to other methods previously presented in the literature, and the same
                                                     identification procedure can be used to identify the fluid-memory models with or without
                                                     having access to the infinite-frequency added mass coefficient. Therefore, it provides an
                                                     extension that puts the two identification problems into the same framework.

                                                     a)                                                                                            b)
                                                                                3D Visualisation of the Wamit file: semisubow.gdf
                                                                                                                                                                                  Convolution Model DoF 33                              x 10 Potential Damping DoF 33

                                                                                                                                                                         150                                                      2.5

                                                                                                                                                                         140                                                       B
                                                                                                                                                                                                                                   Best FD ident, order 8
                                                                                                                                                                                                                       B [Kg/s]

                                                                  10                                                                                                                                                               1
                                                                   0                                                                                                     100                 K (jw) order 8                       0.5
                                                    Z−axis (m)

                                                                                                                                                                          90                                                       0
                                                                 −10                                                                                                       10
                                                                                                                                                                                                                                                                    0    1
                                                                                                                                                                                        Freq. [rad/s]                                               Freq. [rad/s]
                                                                                                                                                                                                                                        x 10
                                                                                                                                                                                                                                             7   Added Mass DoF 33
                                                                 −20                                                                                                     100                                                      7.5

                                                                 −30                                                                                                      50
                                                                                                                                                    Phase K(jw) [deg]

                                                                                                                                                                                                                       A [Kg]

                                                                 −20                                                                                                       0                                                       6
                                                                       0                                                                      40
                                                                                                                                                                                                                                  5.5      A
                                                                           20                                                            20                              −50
                                                                                      40                                            0                                                                                                      Aest FD indet, order 8
                                                                                                60                      −20                                                                                                                Ainf
                                                                                                               −40                                                      −100                                                      4.5
                                                                                                                                                                             −2         −1           0             1                  −2             −1             0    1
                                                                                                                                                                           10         10           10         10                    10             10           10      10
                                                                                                                                Y−axis (m)
                                                                                  X−axis (m)                                                                                            Freq. [rad/s]                                                Freq. [rad/s]

                                                     figure 2: frequency domain identification of marine structure dynamics base on 2D hydrodynamic data.
                                                     (a) shows the geometry of an offshore rig used as an application example (data from
                                            (b) shows the identification results for the frequency response of the
                                                     fluid memory function in heave based on a parametric model of order 8 without using infinite frequency
                                                     added mass as part of the available data. This figure also shows the reconstruction of added mass and
                                                     damping from the identified parametric model.
                                                                                                                             2008 ANNUAL REPORT        21

                                                            A.1.3 Modelling and control of parametric resonance in marine vessels
                                                            Researchers: T. Perez, T. fossen (Norway) and C. Holden (Norway)
                                                            Parametric resonance is a phenomenon where changes in the model parameters can be
                                                            used to describe rapid build-up of oscillations. This phenomenon has been observed in the
                                                            rolling motion of contemporary ship designs with significant bow flare and raised sterns.
                                                            The resonance can be developed when sailing in a head and stern seas with a wave length
                                                            similar to the length of the ship. In these conditions, the hydrodynamic forces that restore
                                                            the up-right equilibrium of the vessel present a time-varying characteristic, which depends
                                                            on the location of the wave crest along the length of the vessel. This is equivalent to a
                                                            mass-spring-damper system where the stiffness of the spring varies with time. This effect
                                                            results in roll parametric resonance due to changes in the restoring forces, which creates
                                                            a rapid development of roll motion that can reach up to 40deg in just a few roll cycles. This
                                                            phenomenon has been responsible for containers being washed overboard in container
                                                            ships and capsizing fishing vessels.
                                                            This project examines the modelling of this phenomenon and the use of controlled roll
                                                            stabilisation devices to reduce the effect on the vessel motion. In 2008, a parametric
                                                            model previously proposed in the literature has been fitted to experimental data of a scale
                                                            model of a 300m containership, and a model of the coupled ship and an actively controlled
                                                            u-tank stabiliser has been developed.

a)                                                          A.1.4 Experiment design and identification of marine vessels dynamics
                                                            Researchers: T. Perez, G.C. Goodwin, J.C. Agüero, E. Herrero (Spain)
                                                                         and T. Armstrong (Austal Ships)
                                                            The mathematical models of marine vessels can be obtained from first principle or
                                                            analytical modelling. The parameters of the model, however, often need to be estimated
                                                            from full scale trials. The ability of an estimation method to produce good estimates of
                                                            the parameters depends on how much information about the dynamics of the system is
                                                            contained in the data, which in turn depends on the experiment performed. Hence, the
                                                            design of experiments is of paramount importance to obtain accurate model parameters
                                                            and to reduce the time to perform experiments. Some of the hydrodynamic phenomena
                                                            involved in ship dynamic response are too complex for analytical modelling.
b)                                                          Therefore, parts of the model structure may also be determined based on
                                                            analysis of experimental data.

                                                Model       This project examines the design of optimal experiments and the application of system
0     50     100    150     200    250    300         350
                                                            identification or experimental modelling for vessels performing different operations:
                                                Estim       dynamic positioning, manoeuvring at low speed and manoeuvring at high speed.

0     50     100    150      200   250    300        350    In 2008, we focused on a two-stage approach for experiment design, in which we first
                      Time [s]
                                                            collect data from step responses, and then based on the information gathered, appropriate
                                                            signals for parametric model identification are designed. This method has been applied

0     50     100    150     200    250    300         350   to the identification of vessels for positioning and slow speed manoeuvring using data
                                                            of full scale trials of a small fishing vessel. In addition, we have proposed the application
                                                            of statistical methods to select hydrodynamic damping model structures that provides

0     50     100    150      200   250    300         350   best use of the information available. These methods have been applied to simplify the
                      Time [s]
                                                            manoeuvring model of a 130m a novel fast ferry trimaran designed and built by Austal
                                                            Ships, Australia.
figure 3: Time-domain identification of a coupled
4DOf manoeuvring model of a highspeed vessel.
(a) shows Austal’s High-Speed Trimaran Design
(Picture courtesy of Austal Ships, Australia).
(b) shows the identification results of a manoeuvring
parametric model with structure selected via stepwise
regression in the time domain. The figure shows the
model and fullscaletrial responses in velocity: surge
(u), sway (v), roll (p), and yaw (r).
22      ARC Centre of Excellence for Complex Dynamic Systems and Control

A.1.5 System identification for rapid model prototyping of ship training
      simulators (Offshore Simulator Centre AS, Norway.)
Researchers: T. Perez and S. Peder-Berge (Norway)
Ship training simulators are used to improve crew efficiency and thus safety of marine
operations. At the core of any virtual-reality simulator lays a mathematical model that
describes the ship dynamic response to control and environmental forces. When a new
vessel is to be incorporated into a simulator, different type of data for the vessel may be
available to be used for system identification to extract a mathematical model for
dynamic response.
In 2008, CDSC evaluated various identification methods for obtaining parametric models
of vessel response based on frequency-domain data computed by hydrodynamic codes,
and made recommendations to the Offshore Simulator Centre AS as to which method to
use for rapid identification of vessel dynamic response.

figure 4: Ship training simulators are used to improve crew efficiency and thus safety of marine operations. At
the core of any virtual-reality simulator lays a mathematical model that describes the ship dynamic response
to control and environmental forces. When a new vessel is to be incorporated into a simulator, different type
of data of the vessel may be available to be used for system identification to extract a mathematical model for
dynamic response.
                                                                                                                 2008 ANNUAL REPORT              23

A.1.6 Adaptive control of gyroscopes
      for roll stabilisation of marine
      vessels (Halcyon International,
                                             reliable computer control systems ensure
                                             adequate operation over an envelope of
                                             sailing conditions.
                                             This project has examined methods for
                                                                                            tool, a designer is able to have a fast initial
                                                                                            assessment of the design before refining
                                                                                            the performance predictions using the
                                                                                            time-domain simulation code GYROSIM.
Researchers: T. Perez and                    determining the size of the gyros to achieve   Also in 2008, an adaptive control strategy
             P. Steinmann (Halcyon)          a desired level of roll reduction and also     was designed to maximize performance
The use of gyroscopic effects of high        for control design to ensure performance       of the gyrostabiliser in changing
speed spinning masses for the roll           in a range of sailing conditions. A time-      environmental conditions. This strategy
stabilisation of marine structures was       domain simulation tool GYROSIM has been        aims at dealing with a fundamental trade-
proposed over 100 years ago. This            developed, which allows a rapid evaluation     off in gyrostabiliser control design in which
approach was very effective, but limited     of the expected performance                    high performance in roll reduction must be
control and large sizes hindered further     of gyrostabilisers.                            balanced with potential constraint violations
developments. In recent years, there has                                                    when a large group of waves arrive at the
                                             In 2008, a frequency-domain simulation
been significant interest in revitalising                                                   vessel. The results have been tested in
                                             tool GYRODESIGN has been developed for
gyrostabilisers due to improvements in                                                      simulation and an experimental prototype
                                             a rapid gyro sizing and to provide design
materials, bearings, and lubricants, which                                                  is being built for further testing in real
                                             information related to mechanical and
have contributed to fast spinning devices                                                   sea conditions.
                                             hydraulic component specification. With this
and size reduction. In addition, fast and

                                                                                            figure 5: Halcyon’s twin gyrostabiliser. A
                                                                                            gyrostabiliser consists of one or more spinning
                                                                                            wheels whose gyroscopic-induced forces are used to
                                                                                            counteract the forces induced by the waves on a ship,
                                                                                            and thus reduce motion. In a roll gyrostabiliser, the
                                                                                            spinning wheel is positioned such that the gyroscopic
                                                                                            effect reduced the motion in roll. The use of twin
                                                                                            wheels rotating in opposite directions eliminates the
                                                                                            gyroscopic effects in other degrees of freedom than
                                                                                            the one intended to be controlled.

                                                                                            figure 6: Numerical study of roll gyrostabilisation of a
                                                                                            navy patrol boat using Haylcyon’s gyro stabilisers. The
                                                                                            top plot shows the roll angle at zero forward speed in
                                                                                            a 3m sea state as a function of the wave period. The
                                                                                            bottom plot shows the expected percentage of roll
                                                                                            angle reduction (RMS).
24    ARC Centre of Excellence for Complex Dynamic Systems and Control

A.1.7 Control design for optimum power extraction in wave energy converters
      (Robotiker-Tecnalia, Spain)
Researchers: T. Perez and M. Santos-Mujica (Spain)
The search for renewable energy resources has revitalised the interest in devices for
wave energy conversion. Wave energy converters (WEC) extract energy from the motion
induced by the waves on particular hull designs. In order to maximize the extracted
energy, the design of a control below to regulate the loads of the power take off element is
of great importance.
In this project, we look at the control aspects of WEC. In particular, we are studying
fundamental performance limitations that can affect potential control system designs
for maximum energy extraction. The preliminary results are exciting and we thus plan
to submit an ARC Discovery Project dedicated to this topic in 2009.
In 2008, we have analysed a particular WEC device, which uses a large gyroscope
mounted on a floating platform. This device extracts energy from precession motion
induced on the gyroscope as a result of the pitching motion of the structure. Based on a
model obtained from a combination of hydrodynamic computations and experimental tests,
we have computed an upper bound on the expected power to be extracted in a various sea
states. We also have guided ROBOTIKER in the development of a time-domain simulation
package, and designed an adaptive precession torque control that aims at optimising the
extracted power.

a)                                                                        b)
                                                                               Power Extraction−Head Seas
                                                                                             Ideal PTO
                                                                                             Gyro Bg = 15, 35, 65 x sqrt(4 Ig Cg)
                                                                                                                                    figure 7: Performance analysis of a wave energy
                                                                                                                                    conversion device with different control designs.
                                                                                                                                    figure (a) shows an schematic of a particular wave
                                                                                                                                    energy converter that uses a power take off element
                                                    Pa/ζ2 [kW/m2]

                                                                                                                                    (PTO) to extract energy from the wave-induced
                                                                                                                                    pitch motion. figure (b) shows the estimated power
                                                                    150                                                             extraction per unit of wave amplitude as a function of
                                                                                                                                    the wave frequency. The ideal PTO is an upper bound
                                                                    100                                                             on performance estimated using only hydrodynamic
                                                                                                                                    characteristics. The use of different control designs
                                                                    50                                                              and a particular PTO was used to evaluate the
                                                                                                                                    potential gains of control adaptation to changes
                                                                     0                                                              in the dominant wave frequency.

A.1.8 Constrained predictive control of ship fin stabilizers to prevent dynamic stall
Researchers: T. Perez and G.C. Goodwin
In moderate to high sea states, the effectiveness of ship fin stabilisers can severely
deteriorate due to nonlinear effects arising from unsteady hydrodynamic characteristics
of the fins: dynamic stall. These nonlinear effects take the form of a hysteresis, and they
become very significant when the effective angle of attack of the fins exceeds a certain
threshold angle. Dynamic stall can result in a complete loss of control action depending
on how much the fins exceed the threshold angle. When this is detected, it is common to
reduce the gain of the controller that commands the fins. This approach is cautious and
tends to reduce performance when the conditions leading to dynamic stall disappear. An
alternative approach for preventing the effects while keeping high performance consists
of estimating the effective angle of attack and set a conservative constraint on it as part
of the control objectives. In our work, we have investigated the latter approach, and here
proposed the use of a model predictive control (MPC) to prevent the development of these
nonlinear effects by considering constraints on both the mechanical angle of the fins and
the effective angle of attack.
                                                                 2008 ANNUAL REPORT         25

Project Leader: J.H. Braslavsky
                   G.J. Adams, J.-C. Agüero G.C. Goodwin, B. Godoy (Student),
                   and A.J. Rojas
External Academic Collaborators:
J.T. Gravdahl (Norwegian University of Science and Technology, Norway),
D. Ugryumova (Student, University of Twente, The Netherlands)
External Industrial Collaborators:
D. Boggs (BHP-Billiton, Perth, Australia)
M. Downey (BHP-Billiton, Newcastle, Australia)
J. Lee (BHP-Billiton, Newcastle, Australia)
A. Maddever (BHP-Billiton, Perth, Australia)
R. Turner (BHP-Billiton, Newcastle, Australia)
This project is funded by a partnership of the Centre with BHP Billiton, and deals with
the development of new technologies using mathematical modelling, and state-of-the-art
model-based control and estimation tools. The project currently encompasses three
n	   Modelling and control of copper heap bioleaching processes,
n	   Sferics reduction in electromagnetic mineral exploration,
n	   Cogeneration at WestVAMP

A.2.1 Modelling and control of copper heap bioleaching processes
Researchers:          J.H. Braslavsky and B. Godoy
This sub-project focuses on the development of mathematical models and control
strategies for heap bioleaching processes for the extraction of copper from sulphide
minerals. Heap bioleaching is of increasing interest in the mining industry to recover
metals from secondary ores. See CDSC Annual Reports 2003-2007 for more
background information.
In previous work, we have suggested the use of feedback control to improve the rate of
mineral extraction based on linearized models around nominal trajectories of the output of
interest. In 2008, this work has been refined in a comparative study between two feedback
approaches: Model Predictive Control (MPC) and Extremum Seeking Control (ESC).
Previously obtained linearized models were used to design an MPC strategy incorporating
input constraints. ESC was tuned, without requiring a process model, to maximise copper
extraction rate using aeration rate.
Simulation results show that similar copper extraction rates can be obtained using either
strategy. The extraction rates improvements with respect to an optimised fixed set-point
strategy are between 4-5%, which correspond to 840-1,050 extra tonnes a year for a
small bioleaching facility. Both feedback strategies improve robustness with respect to
model uncertainties. A paper communicating these results has been presented in the
17th IfAC World Congress in Seoul, South Korea.
This project has been the topic of research of Boris Godoy for his PhD thesis, which
was awarded in October 2008.
26       ARC Centre of Excellence for Complex Dynamic Systems and Control

A    Nominal
         u0             BHPB
                                  Estimated nominal
                                                                Manipulated Process
                                                                    uk                         yk

         Manipulated Process
     +        u                         y + −
                                                         +                              LPF           ×
                                                                     z− 1
              ∆u        Linear         ∆y
                        MPC                                     excitation      A sin(ω k )

figure 8: MPC application to heap bioleaching. A        figure 9: ESC feedback strategy schematics.
high complexity model is used to estimate linearised    Adaptation is applied on the aeration rate setpoint
models around nominal trajectories. MPC is tuned        based on measurements of copper concentration is
to track increments in heap average temperature by      leached solution.
manipulating aeration rate and raffinate influx.
In a real implementation the model could be
retuned adaptively.

figure 10: Comparison of extraction rates improvements for MPC and ESC with respect to an optimised fixed
setpoint strategy. MPC requires less control efforts, while ESC presents slightly faster recovery.
                                                                                                                                                     2008 ANNUAL REPORT           27

  A.2.2 Sferics reduction in
        electromagnetic mineral
  Researchers: J.-C. Agüero,
                                                                             In 2007, broadband noise reduction was
                                                                             achieved using separate models for the
                                                                             low and high frequency ranges. Progress
                                                                             in 2008 has included the development of
                                                                                                                                those obtained using alternative methods,
                                                                                                                                with similar performance on experimental
                                                                                                                                data. The results indicate that the technique
                                                                                                                                can be successfully applied in this case
               J.H. Braslavsky, M. Downey,                                   a single, unbiased model for broadband             and shows potential for other applications.
               G.C. Goodwin, K. Lau,                                         (4 Hz-1 kHz) multinode noise cancellation.         See Lau, Braslavsky, Aguero and Goodwin
               J.B. Lee, A. Maddever,                                        A Matlab implementation of the model               (2008) in Conference Papers.
               P. Turner and                                                 estimation and noise cancellation
                                                                                                                                future work is planned for the extension of
               D. Ugryumova                                                  algorithms has also been written for
                                                                                                                                our noise cancellation techniques to other
                                                                             BHP-Billiton. Other work includes the
  This is a joint project with BHP-Billiton                                                                                     types of magnetic field sensors currently
                                                                             estimation of a model for new sensors.
  Exploration and Mining in Perth                                                                                               being tested by BHP-Billiton.
                                                                             This model can be used for compensation
  (previously in Newcastle).                                                 of the sensor response.
  The aim of this industry project is the                                    In March, Diana Ugryumova, a Masters               A.2.3 Co-generation at WestVAMP
  reduction of sferics noise in mineral                                      Student from the University of Twente in           Researchers: G.J. Adams, G.C. Goodwin,
  exploration using Geoferret, an Australian                                 The Netherlands, visited for six months to                      J.T. Gravdahl (Norway)
  designed and developed electromagnetic                                     complete the practical component of her                         and A.J. Rojas
  exploration system. The exploration                                        course. During her visit she worked on the
  technique employed by Geoferret relies                                                                                        This project is aimed at improvements
                                                                             application of our modelling techniques
  on the induction of currents in the earth                                                                                     to the control for the Westcliff Vent Air
                                                                             to impedance estimation in another
  followed by the measurement of the                                                                                            Methane Project facility near Appin, NSW.
                                                                             electromagnetic exploration method:
  magnetic field generated by the induced                                                                                       The details are subject to a confidentiality
                                                                             magnetotelluric sounding.
  currents. The reduction of sferics noise,                                                                                     agreement, but overall the facility aims to
  electromagnetic noise originating                                          In July, a paper on the application of             generate power from the methane that is
  from lightning storms, is central to the                                   errors-in-variables techniques to sferics          present in air vented from underground
  improvement of signal to noise ratio for                                   attenuation was presented at the IfAC              coal mines. Work in 2008 focussed on
  the detection of deeper ore bodies. See                                    World Congress 2008 in Seoul, Korea.               deriving a model of the process, and
  CDSC Annual Report 2006 for more                                           The technique is used to estimate a model          control improvements are being
  background information.                                                    which is deployed for noise cancellation.          investigated using this model.
                                                                             The estimated model agrees well with

                                      Noise cancellation results − PSD                                                      Noise cancellation results − Coh ZX
                −8                                                                                    0.3

                −9                                                                                    0.2

PSD (V /Hz)


                −11                                                                                      0                  1                  2                 3                    4
              10                                                                                       10                 10                10                 10                10

                                                                                                                                           Coh ZY
              10                                                                                      0.8
                                                                                                      0.6                                                                  residual

                      0           1                   2                  3                4
                   10           10                 10               10               10
                                                   f (Hz)

                                                                                                         0                  1                  2                 3                    4
                                                                                                       10                 10                10                 10                10
                                                                                                                                            f (Hz)

               figure 11: Noise cancellation results. Power spectral density of                         figure 12: Noise cancellation results. Coherence between the Z and X
               the measured sferics noise Z and the residual after performing                           components and the Z and Y components of the sferics noise before and
               noise cancellation. The 50 Hz powerline harmonics have not                               after performing noise cancellation. The coherence has been reduced to a
               been removed.                                                                            negligible level at most frequencies between 4 Hz and 1 kHz. This indicates
                                                                                                        that almost all of the correlated noise has been eliminated.
28    ARC Centre of Excellence for Complex Dynamic Systems and Control

Project Leader: T. Perez
                                                     uses the estimates of one stage to initialise
                                                     the subsequent stage. This algorithm is,
                                                     thus, self-initialised, and the resulting
                                                     estimates also give the option of initialising
Researchers:      G.C. Goodwin, M. fu,               the first proposed algorithm for a further
                  M.M. Zhang,                        refinement of the parameter values. figure
                  B. Godoy (Student),                13 shows a particular realisation of a
                  X. Tai (Student)                   simulated metal price and future contracts.
                                                     figure 14 shows the estimates of the main
External Academic Collaborators:                     parameters of the model corresponding to
K. Barbosa (National Laboratory for                  100 different realisations of simulated data
Scientific Computing, Rio de Janeiro, Brazil)        using the self-initialising algorithm.
External Industrial Collaborators:                   The resulting price models are a critical
P.M. Stone (BHP-Billiton)                            input into BHP Billiton’s models for
M. Menabde (BHP-Billiton)                            understanding the full value of resource
The integrated mine planning project                 assets and for recognizing the value of
aims at developing tools for optimisation            future embedded management options.
of planning and operation of mines. In               The latter are a direct outcome of the
2008, CDSC worked on three particular                level of equilibrium price uncertainty. In
aspects of the problem. The first aspect             an ongoing research robust methods for
is parameter estimation of models for                extending the parameter identification to a
commodity prices. The second aspect                  multi-commodity case are being researched.
concerns the development of planning tools           The second aspect of mine planning looked
and algorithms to optimise both strategic            into in 2008 was mining phase design,
mine plans and infrastructure investment             which is a critical step in long-term mine
decisions. The third aspect is the optimal           planning process and heavily affects the
dispatching scheduling of truck/trains               net present value of a life-of-mine plan. In
within a mine.                                       BHP Billiton’s proprietary mine planning
The major focus of commodity price model             tool, the mine planning process are divided
estimation was to develop algorithms for             into a few sub processes which are
parameter estimation. Commodity price                optimised separately. The mining phases
modelling is normally approached in                  are developed from an “optimal” block
terms of structural time-series models, in           extraction sequence from a prior “block
which the different components (states)              aggregates” scheduling using a fuzzy
have a clear financial interpretation. In            clustering algorithm. The mining phases
particular, we looked at the parameter and           developed in this way sometimes are not
state estimation of two-factor models.               practical, and manual post-processing is
These models use a trend towards an                  required to ensure practical for mining
equilibrium price (long term component)              operation. furthermore, the existing mining
and a reversion to such a trend (short               phase design method does not optimise
term component), which represents the                the net present value.
difference between the current commodity             To address these issues, we have developed
price and the equilibrium price.                     a heuristic automatic method to enforce the
We have developed two different                      practical mining access constraints. We have
algorithms based on the Maximum-                     embedded the method into a framework
likelihood approach. The first algorithm             of meta-heuristic optimisation. Testing has
uses a full parameterisation, and thus               demonstrated that the new approach can
results in a non-linear parameter                    improve the net present value of existing
optimisation problem. This algorithm                 phase designs by the order of 5%.
provides high quality estimates provided             Also some preliminary investigation has
that a good initial guess of the value of            been carried out on the problem of optimal
the parameters is available. The second              train/truck dispatch using approximate
algorithm solves a sequence of estimation            dynamic programming and discrete event
problems that increase in complexity and             dynamic simulation.
                                                                                                                                       2008 ANNUAL REPORT        29

       figure 13: Sample realisation of a commodity
 spot and future contract prices based on simulated
     data. The data is generated using a model with
      parameters estimated from market data of the
                                                                                                               Simulated Data
                           London Metal Exchange.           Spot Price       2.5

                                                                                                                                                   1 month
                                                                               2                                                                   5 months
                                                                                                                                                   9 months
                                                                                   0        50              100          150                     200 months
                                                                                                                                                   13          250
                                                                                                              Time [Weeks]                         17 months
                                                                                                                                                   24 months

                                                            Log Spot Price




                                                                                   0        50              100          150                     200           250
                                                                                                              Time [Weeks]

                                                                                                 Intitial Estimates Usign Spot & Future Prices
     figure 14: Prediction error parameter estimation                        20
results over 100 simulated realizations of commodity
   prices and future contracts. This figure shows the                        18

estimates and true values of the main parameters of                          16
  a 2-factor commodity price model. The parameters
   are the strength of mean reversion (Kappa), trend                         14
  (mu), and the short and long term initial conditions
   of the states of the model. The analysis based on                         12
 simulated data is used to evaluate the properties of                          0       10   20      30       40      50       60       70        80     90     100
    the estimation algorithm; that is, it is used to gain                                                         Realization
              confidence in the results with real data.

                                                                         0             10   20      30       40      50       60       70        80     90     100

                                                                               0       10   20      30       40      50       60       70        80     90     100


                                                                         0             10   20      30       40      50       60       70        80     90     100
          30                      ARC Centre of Excellence for Complex Dynamic Systems and Control

                                                                                 Progress in 2008 occurred in the
                                                                                 following areas:
                                                                                 n	     Tools for economic optimisation,
                                                                                        including data reconciliation strategies
                                                                                                                                                       future work on next-generation model-
                                                                                                                                                       based control tools for CPO will investigate:
                                                                                                                                                       n	     Derivation of MPC tuning parameters
                                                                                                                                                              from parameters of existing (non-MPC)
         Project Leader: G.J. Adams
                                                                                        and interfaces for steady-state                                       control.
         Researchers:                         N. Germyn (Student),                      optimisation of linear and non-linear
                                              G.C. Goodwin, A.M.                                                                                       n	     Robust MPC, based on the algorithm
                                                                                        objectives, were developed and
                                              Medioli, R.H. Middleton,                                                                                        of Løvaas, Seron and Goodwin (see
                                                                                        integrated into CPO (this was completed
                                              M.M. Seron and                                                                                                  Journal Publication)
                                                                                        as an honours student project by Nic
                                              J.S. Welsh                                Germyn).                                                       n	     Testing of the CPOmpc scheduling tool
         External Academic Collaborators:                                                                                                                     applied to non-linear processes.
                                                                                 n	     Code changes in the CPOmpc tool
         D. francois (Université Catholique                                             have resulted in major execution speed                         n	     Extending and integrating the economic
         de Louvain, Belgium)                                                           improvements being attained, and                                      optimisation features as part of a
         External Industrial Collaborators:                                             setup modifications allow on-line model                               “dynamic process optimisation” loop;
         P. farragher (Matrikon)                                                        changes to occur smoothly.                                            some extensions include strategies
         R. Thomas (Matrikon)                                                                                                                                 for extremal seeking, and checks for
                                                                                 n	     Handling of non-linear processes via the
                                                                                                                                                              infeasible setpoints.
                                                                                        appropriate scheduling of multiple linear
         A.4.1 Next generation model based                                              controllers was made possible by the                           n	     A major case study is planned for 2009,
               control tools for CPO                                                    development of an MPC Scheduler.                                      which will bring together the features of
         Researchers:                       G.J. Adams, N. Germyn,                                                                                            non-linear multivariable control, system
                                                                                 n	     Matrikon are developing a CPOmpc
                                            G.C. Goodwin, A.M. Medioli,                                                                                       identification and economic optimisation
                                                                                        control solution platform with an
                                            R.H. Middleton, M.M. Seron,                                                                                       into the control and optimisation of a
                                                                                        overseas company, which will involve
                                            R. Thomas and J.S. Welsh                                                                                          simulated nutating grinding mill.
                                                                                        the control of non-linear systems via
         The aim of this project is to deliver to                                       multiple linear regions.
         Matrikon process control tools that allow                               n	     Decoupling strategies that are intrinsic
         n	             appropriate handling of complex,                                to the internal QP objective function
                        nonlinear and heterogeneous processes;                          of the CPOmpc tool were integrated
                                                                                        and tested. Some results are shown in
         n	             robust and easy-to-use system                                   figure 15 and figure 16.
                        identification; and
         n	             economic optimisation of process

                        1.5                                                                                                         0.4
                                                                                                                                                                                             u1 (original)
     Output amplitude

                                                                                                                                                                                             u1 (decoupled)
                                                                                                                 Input amplitude

                         1                                                                                                          0.2

                        0.5                                                          r1                                              0
                                                                                     y1 (original)
                                                                                     y1 (decoupled)
                         0                                                                                                         −0.2
                              0        50     100     150      200       250   300       350      400                                     0     50      100       150      200       250   300    350     400

                        1.5                                                                                                         0.2
                                                                                                                                                                                             u2 (original)
Output amplitude

                         1                                                                                                          0.1                                                      u2 (decoupled)
                                                                                                                 Input amplitude

                        0.5                                                                                                          0
                         0                                                           y2 (original)                                 −0.1
                                                                                     y2 (decoupled)
                   −0.5                                                                                                            −0.2
                              0        50     100     150      200       250   300       350      400                                     0     50      100       150      200       250   300    350     400
                                                            time (sec)                                                                                                  time (sec)

                              figure 15: Coupled (green) and decoupled (blue) setpoint                                                    figure 16: The associated MV moves in the presence of
                              change responses in the presence of MV constraints.                                                         MV constraints.
                                                                                                                                                     2008 ANNUAL REPORT        31

                                                                A.4.2 Next generation model based control tools for ProcessMORE
                                                                Researchers: P. farragher, D. francois (Belgium) and A.M. Medioli
                                                                This project involves the development of algorithms for automatically generating downtime
                                                                causes from alarm sets. CDSC work in 2008 involved:
                                                                n	                     A thorough investigation of alarm data in relation to physical plant layout and
                                                                n	                     A study of possible or expected classification performance, investigating issues such as
                                                                                       new alarms and downtimes being added over time, irrelevant alarms etc.
                                                                n	                     A visit by Damien francois (Belgium), who developed a simple classifier in 2007.
                                                                n	                     Modification and enhancement of Damien’s classifier, and its application to new data.
                                                                n	                     Moving towards a system that may be used on-line, to adapt to changing alarm sets and
                                                                                       plant configuration.
                                                                The performance of the original “batch” and static classifier and the latest classifier is
                                                                summarised in figure 17, where the “Adaptive Centroid” method would be most useful in a
                                                                real (on-line) system; results show that classification performance for the Adaptive Centroid
                                                                method, when the system suggests up to seven possible causes, has over 50% accuracy.

                                                                                                                                                          Batch Centroid
                                                                                        90                                                                Adaptive Centroid
                                                                                                                                                          Static Centroid
                                                   Correct Classification Percentage







figure 17: Comparison of classifier performance
      for different centroid generation methods.                                         0
                                                                                             1            2              3            4            5              6             7
                                                                                                                             Number of Predictions

                                                                Items to work on in 2009 include:
                                                                n	                     Strategies for culling irrelevant data.
                                                                n	                     Applying the system to other alarm data.
                                                                n	                     Different handling of data pre and post-downtime.
                                                                n	                     Better algorithms to try to improve performance.
32    ARC Centre of Excellence for Complex Dynamic Systems and Control

Project Leader: G.J. Adams
Researchers:      B.J. Burke (Student, CSR Sugar), G.C. Goodwin,
                  A. Rayner (Student), A.J. Rojas and B. Sims
External Academic Collaborator:
J.T. Gravdahl (Norwegian University of Science and Technology, Norway)
External Industrial Collaborator:
R.D. Peirce (CSR Sugar)

A.5.1 Constrained, multi-variable control of an integrated sugar
      mill system for economic enhancement
Researchers: G. Adams, B. Burke, G. Goodwin, J. Gravdahl
             (Norway), R. Peirce, A. Rojas
In today’s carbon – and energy-conscious economy, manufacturers are searching for ways
to maximise returns and minimise wastage. This is certainly true in CSR Sugar, where at
Pioneer Mill (near Ayr, North Queensland) a co-generation plant has been installed. This
plant uses waste cane fibre (bagasse) from a number of mills to create steam for both
sugar milling and for electricity generation. If steam is used for sugar milling efficiently,
CSR can export 50MW of power to the local grid, gaining an income stream from a
waste product.
This project aims to study energy and steam use in sugar processing at Pioneer. The
multi-effect evaporators are core units in the process, where the most efficient evaporation
of water from syrup occurs. Proper control of the sugar content (brix) coming out of the
evaporators, and coordination of steam/energy use with other sections of the mill, are
essential to minimise energy losses.
The current control of the brix suffers from periodic disturbances (90 minute periods),
as well as quicker oscillatory disturbances (10 minute period) caused by addition of
water to the syrup upstream during mill stoppages. CSR and CDSC have studied
evaporator operation, and have determined the cause of the quicker oscillation to be a
type of non-linear flow reduction through the outlet valves. With strategies in place to
reduce this non-linear flow effect in evaporator “5A”, oscillations are reduced from those
seen in figure 18 to those in figure 19. Evaporator “5B” still suffers from these oscillations
in figure 19.
Investigation continues into the cause of the slower (90 minute) oscillations. These
oscillations are most likely to be due to the steam demand from the batch crystallization
pans. Signal analysis is being done on selected process variables to uncover significant
                                                                                                                        2008 ANNUAL REPORT             33

                                                  A.5.2 CSR brake van control
                                                  Researchers: A. Rayner and B. Sims
                                                  Sugar cane trains use a wagon at the end of the train, called a brake van, to control train
                                                  braking. The main aim of the brake van is to keep the couplings between each of the cane
                                                  bins in tension. Once the couplings go into compression, derailments can occur, especially
                                                  when the bins are empty. The operation of the brake van is via radio link from the
                                                  locomotive, and consists of a numbered dial with increasing amounts of brake pressure.
                                                  A park brake can also be independently applied.
                                                  CSR is looking to improve the control of the brake van with an automated process (i.e. the
                                                  brake van automatically selecting the appropriate level of braking for the conditions), as
                                                  well as implementing a new braking unit. In making these improvements, CSR is looking
                                                  to increase the efficiency of the whole process of transporting cane. By automating the
                                                  braking system they will effectively save money on driver training, reduce fuel needs (since
                                                  the brakes are being used more effectively) and cut down on the number of replacements
                                                  to the brake pads. Replacing the current brake type with an electrical braking system will
                                                  save further, since electrical systems brake much faster than the current system.
                                                  Work performed by CDSC on this project so far involves initial modelling of the forces
                                                  involved in the carriages/couplings. further work will concentrate on this modelling aspect,
                                                  and it is envisaged that systems which incorporate GPS information as well may be useful.

figure 18: Examples of poor brix control behaviour in                               figure 19: Improved brix control in 5A (red) compared with 5B (blue).
evaporators 5A (red) and 5B (blue).
34     ARC Centre of Excellence for Complex Dynamic Systems and Control

Project Leader: J.S. Welsh
                                                               A.7 HATCH (INDUSTRIAL
                                                                   CONTROL) (INDUSTRIAL
                                                                                                                                                                                                             A.8 BOEING RESEARCH AND
                                                                                                                                                                                                                 TECHNOLOGY, AUSTRALIA
                                                                                                                                                                                                                 (INDUSTRIAL AffILIATE)
                                                                                                                                                                                                             Project Leader: J.S. Welsh
Researchers:       D. Allingham and
                                                               Project Leader: G.C. Goodwin
                   J.S. Welsh                                                                                                                                                                                Researchers:                                         T. Perez and J.S. Welsh
                                                               Researchers:                                       A. Rojas, C. Renton
External Industrial Collaborators:                                                                                                                                                                           External Industrial Collaborators:
                                                                                                                  (Student), G.C. Goodwin
J. Tusek (Connell Wagner)                                                                                                                                                                                    B.P. Williams (Boeing Research and
                                                               External Industrial Collaborators:                                                                                                            Technology, Australia)
This project, being undertaken in
                                                               T. Domanti (Hatch IAS)                                                                                                                        V. Wheway (Boeing Research and
conjunction with Connell Wagner, is
                                                               G. Wallace (Hatch IAS)                                                                                                                        Technology, Australia)
investigating methods of parameter
estimation for synchronous machines.                           In 2008, the work focussed on cross                                                                                                           This project began in December 2008.
                                                               directional control issues in galvanizing                                                                                                     The details of the project are confidential,
Its aim is to provide validation for
                                                               lines. The project had two streams:                                                                                                           but involve the development of flight
methods already in use, to explore
                                                                                                                                                                                                             systems for autonomous aircraft.
alternative approaches and to investigate                      (i) Rapid estimation of coating thickness
ways of providing error estimates for                          under non-stationary conditions. Here
parameter values.                                              it was found that it was desirable to use
                                                               an estimator having variable memory. (A
Synchronous machines are the primary
                                                               conference paper has been written based
generators of electricity for the power
                                                               on this work).
production industry. Estimation of machine
parameters is a vital field of study, with                     (ii) Combining strip location with double
literature dating back to the 1920s. Many                      sided measurements to give improved
approaches are available, using different                      thickness estimation.
measurement and model regimes. Most
                                                               The latter was principally the work of
current practices are based upon recent
                                                               C. Renton who did this as part of his final
IEEE standards (for example, IEEE Std
                                                               year honours project.
115-1995 and IEEE Std 1110-2002) which
describe in detail both standstill and on-line
tests for synchronous machines.
Our work to date has focussed on standstill                                    ra                     ll                   l fkd1                       l fkd2                         l fd               rfd

frequency response (SSfR) modelling.                                            Id                                                                                                                          Ifd
Here, sinusoidal inputs over a range of                                                                                                                                                                            −1
                                                                                                                                                   l kd1                      l kd 2                        10                                  10                               10                               10
frequencies, approximately from 0.5 mHz to                                                                                                                                                                                                       −1
                                                              Vd                                                                                                                                                                     Vfd
1 kHz, are applied to the machine in a variety                                                                             Lad                                                                                     −2                            −2                               −2                               −2

                                                                                                                                                                                                            10                                  10                               10                               10

of configurations, and the machine response                                                                                                        rkd1                       rkd 2

is measured. Using these responses,                                                                                                                                                                         10
                                                                                                                                                                                                                   −3                            −4
                                                                                                                                                                                                             90                                  90                               90                               90
transfer functions are calculated and the                                                                                                                                                                                                                                                                          60

parameters of an equivalent circuit model                      figure 20                                                                                                                                           60                            60                               60                               30

of the machine, shown in the accompanying                                                                                                                                                                          30                            30                               30                              −30

figure, are then estimated. The measured                                                                                                                                                                                                                                                                          −60
                                                                                                                                                                                                                   0                                 0                                0                           −90
responses along with fits obtained from two                                                                                                                                                                        10
                                                                                                                                                                                                                     −4    −2
                                                                                                                                                                                                                                           10        10
                                                                                                                                                                                                                                                       −4    −2
                                                                                                                                                                                                                                                                            10        10
                                                                                                                                                                                                                                                                                        −4    −2
                                                                                                                                                                                                                                                                                             10      10
                                                                                                                                                                                                                                                                                                         0    2
                                                                                                                                                                                                                                                                                                             10     10
                                                                                                                                                                                                                                                                                                                      −4    −2
                                                                                                                                                                                                                                                                                                                           10      10
                                                                                                                                                                                                                                                                                                                                       0    2

estimation methods are also shown.
                                                                               Zdd                              Zfd                              Zdd0                               pGfd                                       Zff                               Zdf                              Zff0                          pHdf
                                                                   −1                          0                                 −1                                   −1                                           1                             −1                               1                                1
                                                            10                                10                                10                               10                                         10                                  10                               10                               10
The parameters from the SSfR tests will,                                                       −1
                                                                                              10                                                                                                                   0                             −2                               0
                                                                                                                                                                                                            10                                  10                               10
in turn, be used to estimate the machine                                                                                                                                                                                                                                                                           0

                                                                   −2                          −2                                −2                                   −2

                                                            10                                10                                10                               10                                                                                                                                               10

response for on-line step tests, which involve                                                 −3
                                                                                                                                                                                                                   −1                            −3

small changes to the amplitude of the                       10
                                                                   −3                          −4
                                                                                                                                10                               10
                                                                                                                                                                                                                   −2                            −4
                                                             90                                90                                90                               90                                         90                                  90                               90                               90
machine’s power when it is under load (for                                                                                                                            60                                                                                                                                           60

example, when it is connected to the power                         60                          60                                60                                   30                                           60                            60                               60                               30

                                                                                                                                                                       0                                                                                                                                               0
grid). from these tests, time constants are                        30                          30                                30                                  −30                                           30                            30                               30                              −30

estimated which describe the machine’s                                                                                                                               −60                                                                                                                                          −60
                                                                   0                               0                                 0                               −90                                           0                                 0                                0                           −90
response to events such as faults, and how                         10
                                                                     −4    −2
                                                                                         10        10
                                                                                                     −4     −2
                                                                                                                           10        10
                                                                                                                                       −4    −2
                                                                                                                                            10      10
                                                                                                                                                                                           0    2
                                                                                                                                                                                               10                  10
                                                                                                                                                                                                                     −4    −2
                                                                                                                                                                                                                                           10        10
                                                                                                                                                                                                                                                       −4    −2
                                                                                                                                                                                                                                                                            10        10
                                                                                                                                                                                                                                                                                        −4    −2
                                                                                                                                                                                                                                                                                             10      10
                                                                                                                                                                                                                                                                                                         0    2
                                                                                                                                                                                                                                                                                                             10     10
                                                                                                                                                                                                                                                                                                                      −4    −2
                                                                                                                                                                                                                                                                                                                           10      10
                                                                                                                                                                                                                                                                                                                                       0    2

fast it can recover from such events.                                  21:
                                                               figure Zff Plots of transfer function estimates vs. real data for a large synchronous machine.
                                                                                    Zdf           Zff0         pHdf
                                                                   1                           −1                                1                                    1
                                                            10                                10                                10                               10

                                                                   0                           −2                                0
                                                            10                                10                                10

                                                                   −1                          −3                                −1
                                                            10                                10                                10

                                                                   −2                          −4                                −2                                   −1
                                                            10                                10                                10                               10
                                                             90                                90                                90                               90
                                                                   60                          60                                60                                   30

                                                                   30                          30                                30                                  −30

Shared By: