Vol. VIII / 2006 ISSN 1345-8981
ISTECS JOURNAL Science and Technology Policy
[Papers] “Rotational motions, a new observable in seismology” W. Suryanto, H. Igel, A. Cochard, J. Wassermann, U. Schreiber, and A. Velikoseltsev
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 – 13
“Costate Approximation for the State Constrained Optimal Control Problem” Subchan . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 – 24
“Fabrication of Bovine Bone Hydroxyapatite: effect of the material shapes and calcination temperature” M.K. Herliansyah, J.A. Toque, M. Hamdi, A. Ide-Ektessabi, and M.W. Wildan . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 – 33
“PID Parameters Optimization by Using Genetic Algorithm: A Study on Time-delay Systems” A. Mirzal and M. Furukawa
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 – 43
“Enhancement of Agricultural Product Quality Through Management on Controlled Environment” M. A. F. Falah . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 – 53
Published by Institute for Science and Technology Studies (ISTECS) December, 2006
ISTECS JOURNAL VOL. VIII / 2006
1
Editorial Board Chief Editor : Dr. Edi Suharyadi Departemen Fisika UGM, Indonesia & Nagoya University, Japan
Advisory Committee : Dr. Masatoshi Iguchi Bogor Research Station, Indonesia Prof. Yoshiki Mikami Nagaoka University of Technology, Japan Prof. Emer.Dr. Kazutoshi Yamada Chiba Industry Advancement Center, Japan Dr. Mulyanto BATAN, Indonesia Adi Junjungan Mustafa Chiba University, Japan & BAKOSURTANAL, Indonesia
Committee : Dr. Agus Setiyono FKH-IPB, Indonesia Dr. Agus Fanar Syukri P2SMTP-LIPI, Indonesia Dr. Anto Tri Sugiarto KIM-LIPI, Indonesia Dr. Arief Budi Witarto Bioteknologi-LIPI, Indonesia Dr. Budhy Kurniawan FMIPA-UI, Indonesia Dr. Eniya Listiani Dewi P3TM-BPPT, Indonesia Dr. Is Helianti BPPT, Indonesia Dr. Muhammad Budi Setiawan BATAN, Indonesia Dr. Nurul Taufiqu Rochman Fisika-LIPI, Indonesia Dr. Rifaid M. Nur Pengelolaan Limbah Radioaktif, BATAN, Indonesia Dr. Sri Harjanto DTMM, FT-UI, Indonesia Dr. Yopi Biotek-LIPI, Indonesia
ISTECS JOURNAL, Vol. VIII / 2006 -------------------------------------------------------------------------------------------------------------Publisher : Institute for Science and Technology Studies (ISTECS) Jl. Kalibata Selatan No.51 Rt 01/03 Jakarta 12740, Indonesia. Tel / Fax : +62-21-798-4133 http://www.istecs.org, e-mail:
[email protected] -------------------------------------------------------------------------------------------------------------ISTECS JOURNAL VOL. VIII / 2006
2
ISTECS JOURNAL, Vol. VIII (2006) 1 - 13
Rotational motions, a new observable in seismology Suryanto, W.1, H. Igel1, A. Cochard1,2, J. Wassermann1, U. Schreiber3, A. Velikoseltsev3 1
Department of Earth and Environmental Sciences, Geophysics Section, Ludwig-Maximilians Universität München, Theresienstr 41, D-80333, Germany. 2
now at: Ecole et Observatoire des Sciences de la Terre, Institut de Physique du Globe, 5 rue René Descartes, F-Strasbourg Cedex, France 3
Forschungseinrichtung Satellitengeodäsie, Technical University of Munich, Fundamentalstation Wettzell, Sackenriederstr. 25, D-93444 Kötzting, Germany
ABSTRACT The rotational part of earthquake-induced ground motion has basically been ignored in the past decades, compared to the substantial research in observing, processing and inverting translational ground motions. The ring laser has been proven to have consistent rotational result and the capacity to measure the earthquake ranging from several hundred of kilometers distance to some thousand kilometers with a broad range of magnitudes. The consistency is based on the theoretical relation between transverse acceleration and rotation rate: Assuming plane wave propagation they should be in phase and their amplitude ratio proportional to horizontal phase velocity. This suggests that collocated measurements of translations and rotations allow the determination of Love wave dispersion curves, information that is otherwise only accessible through seismic array observations or additional strain measurements.
In this study, we present a novel
approach to estimate frequency-dependent horizontal phase velocities by averaging spectral ratios of translations and rotations from several earthquakes of varying epicentral distances and frequency content. Comparison with theoretical Love-wave dispersion for a spherically symmetric Earth model suggests that the point measurements capture well the expected dispersive behavior. Such point-measurements with additional rotation sensors may proof to be useful for very sparse or single-station networks (e.g. in oceanography or planetary seismology) once appropriate sensors are available. Keywords: ring laser, surface waves, seismic rotations, phase velocity, Love wave dispersion
1. INTRODUCTION The effect of the rotational ground motion has attracted the consideration of only very few seismologists in last decades. Although several classical textbooks (e.g., Davidson, 1927; Gutenberg, 1959; Hobbs, 1907; Imamura, 1937) have mentioned that the effect of rotational motion induced by earthquakes can be observed on the ground surface and some structures, the instrument for measuring ground rotational motion has not yet met the accuracy needed in seismological applications. In earthquake engineering, the rotational ground motion effect has also been recognized for causing structural damage especially for long structures such as bridges and pipelines or transmission systems. In the beginning, those effects were supposed to be due to the asymmetry of the structure or building. Yet, recent studies show that even symmetrical buildings would also be excited into rotational modes Science and Technology Policy ISTECS Journal Vol. VIII/2006
1
about a vertical axis (e.g., Newmark, 1969). Figure 1 shows various effects of rotational motions induced by earthquakes on a tombstone induced by South-central Illinois earthquake November 9, 1968 (Top figures) and M7.0 Miyagi-Oki earthquake of May 26, 2003 (Bottom figures).
Figure 1. Various rotational effect on tombstone induced by earthquake. Top figure: Overturned tombstone after South-central Illinois earthquake November 9, 1968 (Gordonet al. , 1970). Top left: Clockwise rotated tombstone at Campground Cemetery. Top right: Counter-clockwise rotated tombstone at Rector Cemetery. Bottom figures: Rotated tombstone after M7.0 Miyagi-Oki earthquake of May 26, 2003 (Photo Courtesy of The Disaster Control Research Center, Graduate School of Engineering, Tohuku University).
Bouchon & Aki (1982) simulated rotational ground motion near earthquake faults buried in homogeneous layered media for strike-slip and dip-slip fault models. They showed that the maximum rotational velocity produced by a buried 30 km long strike-slip fault with slip of 1 m is 1.5 x 10-3 rad/s. This value is indeed small compared with the amplitude of the translational motion, but in the sense of rotational motion induced by earthquake, those value is very big. They also conclude based on this simulation that rotation, strain and tilt are closely related to ground velocity and the phase velocities associated with strong ground motion are controlled by the rupture velocity and the basement rock shear wave velocity.
Science and Technology Policy ISTECS Journal Vol. VIII/2006
2
Figure 2. Ring laser measurement principle. Two counter-rotating laser beam interfere to generate a beating whien the system rotates with respect to the normal. The beating frequency is directly proportional to rotation rate.
Although rotational sensors were available in the last years (especially for navigation purposes), such instruments with appropriate precision for geophysical application were not available until quite recently. The application of the Sagnac effect for sensing the inertial rotation using optical devices was intensively investigated, since the advent of lasers in sixties. There are two approaches to apply the Sagnac effect for rotational measurements, namely active techniques, as in ring laser gyroscopes, and passive techniques, as in fiber-optic interferometers. The first application of a ring laser gyroscope as a rotational sensor in seismology was reported by Stedman et al. (1997). Ring laser detect the Sagnac beat frequency of two counter-propagating laser beams (Figure 2). The beat frequency δ f is directly proportional to the rotation rate Ω around the surface normal of the ring laser system as given by the Sagnac equation (1) where Ω is the rotation, P is the perimeter of the ring, A the area and λ is the laser wavelength. Furthermore, McLeod et al. (1998) gave a detailed analysis of observations with the ring laser CI, installed in the Cashmere cavern, Christchurch, New Zealand. They reported that the phase of rotation determined by CI is partly consistent with that of a collocated standard seismometer record, during the ML5.3 Kaikoura event on 5 September 1996. Fully consistent rotational motions were recorded by a ring laser gyro installed at the fundamental station Wettzell, Germany (Igel et al., 2005). They showed that the rotational motions induced by the teleseismic event was compatible with collocated recordings of transverse acceleration by a standard seismometer, both in amplitude and phase. In earthquake engineering, observations of rotational components of seismic strong motions may be of interest as this type of motion may contribute to the response of structures to earthquake-induced ground shaking (Li et al., 2001). Rotational motions can be derived from measuring translational motion in several locations (array). Most of rotational/torsional studies of ground motion in earthquake Science and Technology Policy ISTECS Journal Vol. VIII/2006
3
engineering are so far still carried out by indirect measurements. Indirect measurements of rotational motions using a seismo(accelero)meter array have been studied by several investigators (e.g., Bodin et al., 1997; Huang, 2003; Li et al., 2001; Niazi, 1986; Oliveira & Bolt, 1989; Singh et al., 1997; Spudich et al., 1995). However, a comparison of array-derived rotation rate and direct measurement from rotational sensors has never been mentioned in the literature, as no appropriate sensor was available. The full benefits of the determination of rotational motion in seismology are still under investigation. Recent result suggest that the horizontal phase velocity can be estimated by analyzing a collocated recording of rotation and translation ground motion data (Igel et al., 2005). The conventional procedure to estimate the phase velocity is by using array measurements. Using this technique, it would be more cost efficient since it only needs one three component seismometer and a rotational sensor. With array measurements, other properties such as the direction of wave propagation (back azimuth) and the dispersion curve can be derived. If these properties can also be computed from collocated recordings of rotation and translational ground motions, then it may have implication for sparse networks or situations where extremely few or even single-station observations are taken (e.g. in remote areas or planetary seismology). In this paper we will show a specific application of collocated measurements of translations and rotation to determine the Love-wave dispersion curve. The determination of frequency-dependent surface-wave phase velocities has for a long time been one of the most important tools to determine 3-D seismic velocity structure on regional and global scales (e.g., Nataf et al., 1984, 1986; Snieder, 1988ab, and many others). On small scales, near-surface low-velocity structures – crucial for the estimation of hazard-relevant site effects – can be determined using ambient noise measurements (e.g., Milana et al., 1996, Kind et al., 2005). Recently, it was shown that Rayleigh-wave dispersion curves can be derived by correlating long time series of ambient noise (micro-seismicity) and that the velocity structure thus derived can be used to image 3-D structures (Campillo and Paul, 2003; Shapiro and Campillo, 2004; Shapiro et al., 2005). The aforementioned techniques require observations from seismic arrays to recover frequency-dependent propagation times (and thus phase velocities) across the array in the direction of propagation. Standard seismic observations are restricted to three components of translations, despite the fact that the recovery of the complete motion requires the observation of three additional components of rotations and six components of strain (e.g., Aki and Richards, 2002; Trifunac and Todorovska, 2001). In the past years, rotation sensor technology has been improving in a way that may allow the development of routine sensors for three additional rotational motion components useful for seismological purposes (e.g., Schreiber et al., 2005, 2006). Recent observations of local, regional and global wavefields using ring laser technology showed that the rotational measurements are fully consistent with collocated observations of translations (e.g., Igel et al., 2005; Cochard et al., 2006; Igel et al., 2006) following earlier observations of earthquake-induced rotational motions (McLeod et al., 1998; Pancha et al., 2000). Further confirmation of accurate measurements of the new observational component using ring laser technology came through comparison with rotational motions derived from seismic array data (Suryanto et al., 2006) using a classical approach (e.g., Spudich et al. (1995)). A temporary array was installed Science and Technology Policy ISTECS Journal Vol. VIII/2006
4
around the ring laser instrument and direct and array-derived rotations compared for an event with high signal-to-noise ratio. The high correlation-coefficient (0.93 %) and almost identical amplitudes for the two independent rotation measurements observed with entirely different physical principles further indicate that the ring laser indeed measures the rotational motions accurately in a wide frequency range. A simple relationship between transverse acceleration and rotation rate (around a vertical axis) shows that both signals should be in phase and their ratio proportional to horizontal phase velocity. Igel et al. (2005) and Cochard et al. (2006) exploited this relationship to estimate horizontal phase velocities in sliding time windows along the observed time series. Comparison with synthetic traces (rotations and translations) and phase velocities determined in the same way showed good agreement with the observations. These initial results suggested that the determination of Love-wave dispersion curves (and thus information on local 1-D shear velocity structure) may be possible. It is worth noting that a similar relationship between strain and displacements can be used to determine horizontal phase velocities (e.g., Mikumo and Aki, (1964), Gomberg and Agnew (1996)). In this study we present a novel method for the determination of Love-wave phase velocities based on collocated measurements of translations (standard broadband seismometer) and rotations around a vertical axis (observed by a ring laser). Instead of determining phase velocities in the time domain (Igel et al., 2005; Cochard et al., 2006), we average spectral ratios of several earthquakes which allows us to directly determine frequency dependent Love-wave phase velocities and compare them with theoretical predictions for spherically symmetric Earth models. The results are supported by applying the same processing steps to complete 3-D synthetic seismograms for some of the observed events. Table 1. Parameters of the earthquakes used in this study.
Date
Time (UTC)
Lon.
Lat.
Magnitude
Location
21/05/03
18:44:19
003.71E
36.90N
6.8
Northern Algeria
26/05/03
09:24:32
141.45E
38.90N
7.0
Honshu, Japan
06/07/03
19:10:33
026.07E
40.34N
5.7
Turkey
14/08/03
05:14:55
020.74 E
39:19N
6.3
Greece
25/09/03
19:50:06
143.90E
41.77N
8.3
Hokkaido, Japan
27/09/03
18:52:53
087.69E
50.06N
6.6
Southern Siberia
01/10/03
01:03:25
087.68E
50.22N
6.7
Southern Siberia
05/02/04
21:05:24
135.49E
03.58S
7.0
Irian Jaya, Indonesia
27/02/04
02:27:46
003.96W
35.23N
6.3
Gibraltar
2. OBSERVATIONS AND DATA PROCESSING We use translation data from station Wettzell (WET) of the German Regional Seismic Network (GRSN) located in Southern Germany (12o52’44”E, 49o08’39”N). The station is equipped with an STS-2 broadband instrument with a flat response to ground velocity from 8.33 mHz (120 s) to 50 Hz. Science and Technology Policy ISTECS Journal Vol. VIII/2006
5
The data with a sampling rate of 80 Hz are corrected for instrument response, rotated into a local radial-transverse system, and differentiated to obtain transverse acceleration. The rotational data are measured by a ring laser instrument, called “G”, consisting of a He-Ne gas laser with an ultrahigh vacuum quality cavity enclosing an area of 16 m2. The vertical component of rotation rate is recorded by this instrument with a sampling rate of 4 Hz. The instrumental sensitivity of ring lasers is limited by the scale factor and quantum noise processes. For the G ring laser rotation rates as small as 10-10 rad/s/◊Hz can be observed (Schreiber et al., 2003). Further information on the ring laser instrument is given in Schreiber et al. (2005). The ring laser is mounted horizontally in the Geodetic Fundamentalstation Wettzell (about 250 m from the STS-2 seismometer). Given the frequency range (i.e., spatial wavelengths) considered below, we treat the two observations (rotations and translations) as collocated. From a growing event database with translations and rotations (see Igel et al., 2006) we use several regional and global earthquakes in 2003 and 2004 with M > 5.7, listed in Table 1. 3. EXAMPLES OF PHASE DETERMINATION IN THE TIME-DOMAIN We first illustrate the possibility of deriving phase velocities using the time-domain approach pursued by Igel et al. (2005, 2006). In Figure 3a+b, time series of transverse acceleration (gray) and rotation rate (black) are shown for two events, the M6.3, Greece, 14 August, 2003, and the M6.7, Siberia, 1 October, 2003, respectively. The almost identical waveform fit between rotations and translations in both cases illustrate that the assumption of plane wave propagation is appropriate and that information on the horizontal phase velocity should be contained in the ratio between transverse acceleration and rotation rate. An appropriate measure of the fit between two presumably synchronous signals is the zero-lag normalized cross-correlation coefficient. We quantify the time-dependent similarity between rotation rate and transverse acceleration by sliding a time-window (10 s) along the time series and calculate the cross-correlation coefficient that is defined between 0 (no similarity) and 1 (perfect match). If the quality of the waveform fit in a given time window is above a threshold (0.95) we estimate a horizontal phase velocity for this time window by finding the best-fitting velocity in a least-squares sense, as well as the associated variance. These phase velocities and the associated uncertainties are shown for two particular earthquakes in the bottom plots of Figures 3a and 3b for time windows containing the fundamental Love-waves mode. In both cases, the estimated phase velocities are within the expected range of fundamental mode Love-wave phase velocities for spherically symmetric Earth models (3-5 km/s). However, the time-domain representation makes it difficult to extract the frequency dependent behavior of Love waves. Therefore, we introduce an approach in which the phase velocities are directly estimated in the frequency domain.
Science and Technology Policy ISTECS Journal Vol. VIII/2006
6
(a)
(b)
Figure 3. a: Upper trace: Transverse acceleration (gray, left axis) and rotation rate about the vertical axis (black, right axis) for the Greece event, M6.3, 14 August 2003. Bottom trace: Best-fitting horizontal phase velocities as a function of time in a 10 s sliding window. b: Same for the Siberia event, M6.7, 1 October 2003.
4. LOVE WAVE DISPERSION The results above and those reported by Igel et al. (2006) and Cochard et al. (2006) indeed suggest that it should be possible to determine the phase velocities as a function of frequency (dispersion) by calculating the spectral ratios of transverse acceleration and rotation rate for time windows containing the Love-wave trains. For this purpose, the rotation rate is interpolated to the same sampling points as the transverse acceleration and the Love wave train time window isolated and attenuated at the edges with a Gaussian function. Both time series are transformed into the Fourier domain and the ratio of the spectra of rotations and transverse acceleration leads to the Science and Technology Policy ISTECS Journal Vol. VIII/2006
7
frequency-dependent phase velocities aT (ω) = −2c(ω). Ω(ω)
(2)
Because of the oscillatory nature of the individual spectra and spectral ratios we average the ratios from several events assuming that the resulting phase velocities are representative of the same subsurface volume. In addition, we smooth the ratios along the frequency axis using a Savitzky-Golay filter, a low pass filter also known as least square smoothing filter or DISPO (digital smoothing polynomial). The filter is defined as a weighted moving average with weights given as a polynomial of a certain degree; in this case we use degree two (Press et al., 2002).
Figure 4. Frequency-dependent Love wave phase velocities from spectral ratios of synthetic seismograms calculated for the Gibraltar, Papua and Hokkaido events (see text for details) shown as mean values with variances as Gaussian uncertainties. The dashed lines indicate the theoretical fundamental- (lowest dashed line) and higher-mode Love-waves phase velocities (upper dashed lines) obtained for the ak135 Earth model.
We first test the methodology presented above on complete synthetic seismograms (rotations and translations) calculated for one regional (Gibraltar) and two global (Hokkaido and Papua) events using the spectral-element method (Komatitsch and Tromp, 2002a,b) employing a recent 3-D tomographic model (Ritsema and Van Heijst, 2000) and the crust model by Bassin et al. (2000). The sources are modeled as point shear dislocations with source properties from the Harvard Catalogue [http://www.seismology.harvard.edu/]. The resulting spectral ratios were averaged and processed as described above. The frequency-dependent phase velocities are shown in Figure 4 as Gaussians with mean value and variance for each period. We superimpose theoretical predictions of Love-wave dispersion curves for the spherically symmetric ak135 Earth model (Kennett et al., 1995) for the fundamental and the first three higher-order modes. Despite the small number of events the estimated Love-wave phase velocities seem to capture well those predicted for the fundamental modes in a Science and Technology Policy ISTECS Journal Vol. VIII/2006
8
spherically symmetric Earth model. The uncertainties decrease with increasing period. In the frequency (period) window considered, the largest deviation from the predicted values are 6.1% (at period 20s).
Figure 5. Love waves phase velocities derived from observed data shown as mean values (black) and associated uncertainties using grey shading (see text for details) for all event listed in Table 1. The dashed lines indicate fundamental and higher-mode Love-waves phase velocities obtained for the ak135 Earth model.
We calculate stacked spectral ratios for the regional (Greece, Turkey, Gibraltar, Algeria) and global (all other) events listed in Table 1. The results are presented in Figure 5 in the same way as the synthetic data shown in Figure 4. Except in the period range around 10 s, the Love-wave phase velocities determined by the spectral ratios are close to the predicted values. The increase in phase velocities towards longer periods is well captured by the stacked global events, with almost constant uncertainty with period. The maximum deviations of the mean values are about 5.8% (at period 100 s). The reasons for any discrepancies may be 3-D heterogeneity, anisotropy, non-planar wavefronts, deviations from the great circle paths, or uncertainties in the observations of translations and rotations (e.g., site effects). Fully understand these potential sources of discrepancies will require a larger database, comparison with array-derived Love-wave dispersion curves and systematic synthetic studies. 5. CONCLUSIONS The aim of this study was to present a novel methodology to derive Love-wave dispersion curves with point measurements of rotations (around a vertical axis) and translations. Frequency-dependent phase velocities are estimated by calculating the ratio between the spectra of transverse acceleration and rotation rate by stacking the ratios of several events. This approach was applied to 3-D synthetic data Science and Technology Policy ISTECS Journal Vol. VIII/2006
9
sets and several regional and global events observed by the collocated ring laser instrument measuring the rotation rate around a vertical axis and a standard broadband sensor located at Wettzell, SE-Germany. Both synthetic and observed dispersion curves match well those predicted for the fundamental mode Love-waves. This indicates that plane-wave theory is appropriate and that the assumption of fundamental mode Love-wave propagation is approximately fulfilled, or that the energy of higher-mode Love waves in the time-windows considered is low. The purpose of this study was primarily to illustrate the concept and show a first application to real observations. Love-wave dispersion curves can be used to derive local 1-D velocity structure and are therefore an important intermediate result for tomographic inversions. Whether the accuracy of the dispersion curves derived with the approach presented here is enough for tomographic purposes remains to be evaluated. We intend to investigate these issues by systematic synthetic studies and analysis of a larger event database. Nevertheless, the results shown here indicate that through additional measurements of accurate rotational signals, wavefield information is accessible that otherwise requires seismic array data. This may be of use when arrays are very sparse or consist of only one station (e.g., oceanography or planetary seismology). However, to make this methodology practically useful for seismology will require the development of an appropriate high-resolution six-component broadband sensor. Efforts are underway to coordinate such developments on an international scale (Evans et al., 2006). ACKNOWLEDGEMENTS This work is supported by the German Ministry of Research and Education (BMBF-Geotechnologien) and German Research Foundation. W.S. was supported by the German Academic Exchange Service (IQN-Georisk), and by the Graduate program THESIS funded by the Bavarian Government. We are grateful to J. Tromp and D. Komatitsch for providing their Spectral Element Method code. The theoretical Love wave dispersion curve was calculated using the program lmtor provided by Wolfgang Friederich (available through the SPICE library, www.spice-rtn.org). REFERENCES [1] Aki, K. and P. Richards, 2002. Quantitative Seismology, 2nd Edition, University Science Books. [2] Bassin, C., G. Laske, and G. Masters, 2000. The Current Limits of Resolution for Surface Wave Tomography in North America, EOS Transactions, AGU 81, F897. [3] Bouchon, M. & Aki, K., 1982. Strain, tilt, and rotation associated with strong ground motion in the vicinity of earthquake faults. Bull. Seism. Soc. Am., 72, 1717-1738. [4] Campillo, M., and A. Paul, 2003. Long range correlations in the diffuse seismic coda, Science, 299, 547-549. [5] Cochard, A., H. Igel, A. Flaws, B. Schuberth, J. Wassermann, and W. Suryanto, 2006. Rotational motions in seismology: theory, observation, simulation, in Teisseyre et al., editor, Earthquake source asymmetry, structural media and rotation effects. Springer Verlag, 391-412. [6] Davidson, C., 1927. The Founders of Seismology. Cambridge university press. Science and Technology Policy ISTECS Journal Vol. VIII/2006
10
[7] Evans J.R., A. Cochard, V. Graizer, B.-S. Huang, K. W. Hudnut, C. R. Hutt, H. Igel, W.H.K. Lee, [8] C.-C. Liu, R. Nigbor, E. Safak, W.U. Savage, U. Schreiber, M. Trifunac, J. Wassermann, and C.-F. Wu, 2006. A workshop on rotational ground motion, Seismol. Res. Lett., submitted. [9] Gomberg, J., and Agnew, D., 1996. The accuracy of seismic estimates of dynamic strains: An evalutation using strainmeter and seismometer data from Pinon flat Observatory, California, Bull. Seismol. Soc. Am., 86, pp. 212-220. [10] Gutenberg, B. (1959). Physics of the Earth's Interior. Academic Press. [11] Hart, G.C., Dijulio, R.M. & Lew, M., 1975. Torsional response of high-rise buildings. ASCE, Journal of the Structural Division, 101, 397-415. [12] Hobbs, W.H., 1907. Earthquakes. Appleton. [13] Igel, H., Cochard., A., Wassermann, J., Flaws., A., Schreiber, U., Velikoseltsev, A., Nguyen, P. D., 2006. Broadband Observations of Earthquake Induced Rotational Ground Motions, Geophys. J. Int., (in press). [14] Igel, H., U. Schreiber, A. Flaws, B. Schuberth, A. Velikoseltsev, and A. Cochard, 2005. Rotational motions induced by the M 8.1 Tokachi-Oki earthquake, September 25, 2003, Geophys. Res. Lett., 32, L08309. [15] Immamura, A., 1937. Theoretical and Applied Seismology. Maruzen. [16] Kennett, B.L.N. Engdahl, E.R. and Buland R., 1995. Constraints on seismic velocities in the Earth from travel times, Geophys J. Int., 122, 108-124 [18] Kind, F., Fäh, D., and Giardini, D., 2005. Array measurements of S-wave velocities from ambient vibrations, Geophys. J. Int, 160, pp. 114-126. [19] Komatitsch, D., and J. Tromp, 2002a. Spectral-element simulations of global seismic wave propagation, Part I: Validation. Geophys. J. Int., 149, 390-412. [20] Komatitsch, D., and J. Tromp, 2002b. Spectral-element simulations of global seismic wave propagation, Part II: 3-D models, oceans, rotation, and gravity. Geophys. J. Int., 150, 303-318. [21] McLeod, D. P., G. E Stedman, T. H. Webb and U. Schreiber, 1998. Comparison of standard and ring laser rotational seismograms, Bull. Seismol. Soc. Am., 88, 1495-1503. [22] Mikumo, T., and K. Aki, 1964. Determination of local phase velocity by intercomparison of seismograms from strain and pendulum instruments, J. Geophys. Res., 69, 721-731. [23] Milana G., Barba, S., del Pezzo, E., Zambonelli, E., 1996. Site response from ambient noise measurements: new perspectives from an array study in Central Italy, Bull. Seis. Soc. Amer., 86, 320-328. [24] Nataf, H.-C., I. Nakanishi, and D.L. Anderson, 1984. Anisotropy and shear velocity heterogeneities in the upper mantle, Geophys. Res. Lett., 11, 109-112. [25] Nataf, H.-C., Nakanishi, and D.L. Anderson, 1986. Measurement of mantle wave velocities and inversion for lateral heterogeneity and anisotropy, III. Inversion, J. Geophys. Res., 91, 7261-7307. [26] Newmark, N.M., 1969. Torsion in symmetrical buildings. Proc. Fourth World Conference on Earthquake Engineering, Santiago, Chile, 3, 19-32. [27] Pancha, A., T. H. Webb, G.E. Stedman, P. McLeod and U. Schreiber, 2000. Ring laser detection of Science and Technology Policy ISTECS Journal Vol. VIII/2006
11
rotations from teleseismic waves, Geophys. Res. Lett., 27, 3553-3556. [28] Press, W. H., Teukolsky, S. A., Vetterling, W. T., and Flannery, B. P., 2002. Numerical Recipes, Cambridge University Press, pp. 650-655. [29] Ritsema, J., and H. J. Van Heijst, 2000. Seismic imaging of structural heterogeneity in Earth’s mantle: Evidence for large-scale mantle flow, Science Progress, 83, 243-259. [30] Schreiber, K. U., A. Velikoseltsev, G. E. Stedman, R. B. Hurst, and T. Klügel, 2003. New applications of very large ring lasers, in Symposium Gyro Technology, edited by H. Sorg, 80– 87. [31] Schreiber, K.U., H. Igel, A. Cochard, A. Velikoseltsev, A. Flaws, B. Schuberth, W. Drewitz, and F. Müller, 2005. The GEOsensor project: rotations - a new observable for seismology. In Rummel et al., editor, Geotechnologies, Springer Verlag. [32] Schreiber, U., G.E. Stedman, H. Igel, and A. Flaws, 2006. Ring laser gyroscopes as rotation sensors for seismic wave studies. In “Earthquake source asymmetry, structural media and rotation effects” eds. Teisseyre et al., Springer Verlag, 377-390. [33] Shapiro, N. M., and Campillo, M., 2004. Emergence of broadband Rayleigh waves from
[34] [35] [36] [37]
correlations of the ambient seismic noise, Geophys. Res. Lett., 31, L07614, doi:10.1029/2004GL019491. Shapiro, N.M., M. Campillo, L. Stehly, M.H. Ritzwoller, 2005. High-resolution surface wave tomography from ambient seismic noise, Science, 307, 1615-1618. Snieder, R., 1988a. Large scale waveform inversion of surface waves for lateral heterogeneity, 1, Theory and numerical examples, J. Geophys. Res., 93, 12055-12065. Snieder, R., 1988b. Large scale waveform inversion of surface waves for lateral heterogeneity, 2, application to surface waves in Europe and the Mediterranean, J. Geophys. Res., 93, 12067-12080. Spudich, P., L. K. Steck, M. Hellweg, J. B. Fletcher, and L. M. Baker, 1995. Transient stresses at Parkfield, California, produced by the M7.4 Landers earthquake of June 28, 1992: Observations from the UPSAR dense seismograph array. J. Geophys. Res., 100, 675-690.
[38] Stedman, G.E. (1997). Ring laser tests of fundamental physics and geophysics. Rev. Progr. Phys., 60, 615Œ688. [39] Suryanto, W., Igel, H., Wassermann, J., Cochard, A., Vollmer, D., Scherbaum, F., Schuberth, B., Schreiber, U., Velikoseltsev, A., 2006. First comparison of array-derived rotational ground motions with direct ring laser measurements, BSSA (in press). [40] Trifunac, M.D., and M. I. Todorovska (2001). A note on the useable dynamic range of accelerographs recording translation. Soil Dyn. and Earth. Eng., 21, 275-286.
Science and Technology Policy ISTECS Journal Vol. VIII/2006
12
Profile Wiwit Suryanto: Geophysic Lab. Staff, Gadjah Mada University, Curretly doing PhD research at the University of Munich, Germany. Heiner Igel: Professor of Geophysic at the University of Munich Germany. Alain Cochard: Researcher at
Ecole et Observatoire des Sciences de la Terre, Institut de Physique du Globe, 5 rue
René Descartes, F-Strasbourg Cedex, France. Joachim Wassermann: Head of the Geophysical Observatory Fuerstenfeldbruck, University of Munich, Germany. Ulrich Schreiber: Professor at the Technical University of Munich. Alex Velikoseltsev: Researcher at the Fundamentalstation Wettzell, Germany.
Science and Technology Policy ISTECS Journal Vol. VIII/2006
13
ISTECS JOURNAL, VIII (2006) 14 – 24
Costate Approximation for the State Constrained Optimal Control Problem Subchan1,2 1
Department of
Mathematics, Institut Teknologi 10 Nopember Surabaya
Kampus ITS Keputih, Sukolilo, Surabaya, Indonesia 2
Department of Aerospace Power and Sensors, Cranfield University Defence College of Management and Technology Defence Academy of the United Kingdom Shrivenham, Swindon SN6 8LA, United Kingdom
ABSTRACT This paper presents an analysis and numerical solutions for costate approximation of the optimal control problem with pure state constraint. The method is based on the direct and indirect methods. The direct approach replaces the continuous time interval with a grid of discrete points, thus approximating it with a finite-dimensional problem, albeit of high dimension. The indirect approach preserves the infinite-dimensional character of the task and uses the theory of optimal control to solve it. The direct method solutions based on the Jacobson's formulation are compared with the indirect method solutions which are based on the Bryson's formulation. Examples are presented that show the inaccurate results of the Jacobson's formulation. These results provide useful insight to the optimal control problem with active pure state constraint. Keywords: Optimal Control; Direct Method; Indirect Method; Pure State Constraint.
1. INTRODUCTION The main concern in the solving optimal control problem by an indirect method approach is defining the initial guesses for the costate variables especially when the state constraint is active. Von Stryk [7] shows that the costate variable can be estimated by the necessary conditions of the discretised problem of the optimal control. He developed the DIRCOL package [8] based on a direct collocation method and it has been used for solving several real-life problems (see von Stryk and Bulirsch [9] and von Stryk and Schlemmer [10]. Grimm and Markl [4] estimated the costate variables using direct multiple shooting methods. Their costate approximation is accurate for the unconstrained problem while it does not work well for the constrained problem. Fahroo and Ross [3] proposed a Legendre pseudospectral method to estimate the costate variables and presented an accurate estimator for the unconstrained problem. Benson [1] proposed a Gauss pseudospectral transcription to solve optimal control problem and used it to approximate costate variables. Again the costate approximation does not give good initial guess for the pure state constrained problem. Science and Technology Policy ISTECS Journal Vol. VIII/2006
14
This paper presents the direct and indirect approaches of computing the costate variables. The Bryson's and Jacobson's formulation are compared and then the DIRCOL package is used to solve the both formulations. 2. OPTIMAL CONTROL PROBLEM The problem is to find an admissible control u (t ) , which minimises the performance index: tf
min J = Φ[ x(t f ), t f ] + ∫ L[ x(t ), u (t ), t ]dt
(1)
t0
u∈U
with respect to the state vector functions:
{
}
Χ = x : [0, t f ] → ℜ n xi , i = 1,K, n, piecewise continuously differentiable ,
and the control vector functions:
{
}
U = u : [0, t f ] → U ⊂ ℜ n u i , i = 1,K, n, piecewise continuous ,
(2)
(3)
subject to the following constraints: x& = f( x(t ), u (t ))
f : ℜ n+ m → ℜ n
(4)
x(0) = x0 ∈ ℜ n
x0 known
(5)
ψ x(t f ), t f = 0 ∈ ℜ p
ψ : ℜ n xℜ + → ℜ p , p ≤ n, t f unknown
(6)
C( x(t ), u (t ) ) ≤ 0 ∈ ℜ q
C : ℜ n+ m → ℜ q
(7)
S(x(t ) ) ≤ 0 ∈ ℜ s
S : ℜn → ℜs
(8)
(
)
The performance index describes a quantitative measure of the performance of the system over time. In aerospace problems, a typical performance index gives an appropriate measure of the quantities such as minimum fuel/energy, optimal time etc. Here Φ : ℜ n+1 → ℜ1 and L : ℜ n + m → ℜ1 are assumed to be sufficiently often continuously differentiable in all arguments. Minimising J with respect to the control function u must be accomplished in a way consistent with the dynamics of the system, whose performance is optimised. In other words, equation (4) is the first fundamental equality constraint. The optimal control u ∗ , when substituted to (4), will produce the optimal state x ∗ , while minimising J . The optimal state x ∗ is further constrained by the boundary conditions (5) and (6), path constraints (7) and (8). Finally, it should be mentioned that (7) or (8) may be active on a subinterval of
(t 0 , t f ) or just at a point. In the former case, the constrained (active) subarc will be characterised by the
Science and Technology Policy ISTECS Journal Vol. VIII/2006
15
entry time t1 and the exit time t 2 with t 0 ≤ t 1 ≤ t 2 ≤ t f .
T L
« V
x D
mg
Figure 1: Definition of missile axes and angles
The functions appearing in (1)-(8) are assumed to be sufficiently continuously differentiable with respect to their arguments. Note that the definition of U allows discontinuities in controls and thus implies corners (cusps) in the states, so that X comprises piecewise smooth functions. This is a practical necessity, as many real-world applications of optimal control involve bang-bang type inputs. Problem (1)-(8) is infinite-dimensional: its solution is not a finite vector of numbers, but a function. For a real-life application it is impossible to guess the optimal function, so a recourse to approximate methods is necessary. They attempt to find a finite-dimensional representation of the solution which is accurate at the nodes of the representation, has acceptable error between the nodes and converges to the true function as the number of nodes tends to infinity, if second order sufficient conditions hold. There are two main approaches to the solution of the problem. The direct approach replaces the continuous time interval with a grid of discrete points, thus approximating it with a finite-dimensional problem, albeit of high dimension (hundreds of discretised variables). The indirect approach preserves the infinite-dimensional character of the task and uses the theory of optimal control to solve it. 3. PROBLEM FORMULATION Trajectory shaping of a missile is an advanced approach to missile guidance which aims at computing the whole trajectory in an optimal way. The mission is to hit a fixed target while minimising the missile exposure to anti-air defences. The problem is to find the trajectory of a generic cruise missile from the assigned initial state to a final state with the minimum altitude along the trajectory. The objective can be formulated by introducing the performance criterion: Science and Technology Policy ISTECS Journal Vol. VIII/2006
16
tf
J = ∫ h dt
(9)
t0
The performance criterion is subject to the equations of motion, which may be written as:
γ& =
T −D L g cos γ sin α + cos α − mV mV V
(10a)
T −D L V& = cos α − sin − g sin γ m m x& = V cos γ
(10b) (10c)
h& = V sin γ where
t
is
the
(10d)
actual
time,
t0 ≤ t ≤ t f
with
t 0 as
the
initial
time
and
tf
as the final time. The state variables are the flight path angle γ , speed V, horizontal position x and altitude h of the missile. The thrust magnitude T and the angle of attack α are the two control variables (see Figure 1). The aerodynamic forces D and L are functions of the altitude h, velocity V and angle of attack α . The following relationships have been assumed: Drag. The drag D is written in the form
D(h,V , α ) =
1 C d ρV 2 S ref 2
(11)
C d = A1α 2 + A2α + A3
(12)
Lift. The lift L is written in the form
1 Cl ρV 2 S ref 2 Cl = B1α + B2
L(h,V , α ) =
(13) (14)
Where ρ is air density given by
ρ = C1h 2 + C 2 h + C3
(15)
and S ref is the reference area of the missile; m denotes the mass and g the gravitational constant. Boundary conditions. The initial and final conditions for the four state variables are specified:
Science and Technology Policy ISTECS Journal Vol. VIII/2006
γ ( 0) = γ 0
γ (t f ) = γ t f
V (0) = V0
V (t f ) = Vt f
x ( 0) = x 0
x(t f ) = xt f
h(0) = h0
h(t f ) = ht f
(16)
17
Table 1: Boundary data, constraint data and physical constants Quantity
Value
Unit
Quantity
Value
Unit
Vmin
200
m/s
Lmax
4
Vmax
310
m/s
A1
-1.9431
M
1005
kg
A2
-0.1499
G
9.81
m/s2
A3
0.2359
Sref
0.3376
m2
B1
21.9
Tmin
1000
N
B2
0
Tmax
6000
N
C1
3.312 10-9
kg m-5
hmin
30
m
C2
-1.142 10-4
kg m-4
Lmin
-4
C3
1.224
kg m-3
Figure 2: Dircol computational results for the state variables of Bryson’s formulation.
In addition, constraints are defined as follows:
•
State path constraints
Science and Technology Policy ISTECS Journal Vol. VIII/2006
Vmin ≤ V ≤ Vmax
(17)
hmin ≤ h
(18) 18
•
Control path constraint Tmin ≤ T ≤ Tmax
•
Mixed state and control constraint L Lmin ≤ ≤ Lmax mg
(19)
(20)
where Lmin and Lmax are normalized, see Table 1.
3. COSTATE APPROXIMATION This section presents an example of the state constrained optimal control problem. In this example, Bryson's and Jacobson's formulation are compared and then the DIRCOL package is used to approximate the costate variables. Jacobson et al. [9] presented a direct adjoining of pure state constraint (18) to the Hamiltonian while Bryson et al. [10] proposed an indirect adjoining of pure state constraints to the Hamiltonian. In Bryson's approach the pure state constraint is differentiated until appears explicitly and then the resulting equation is adjoined to the Hamiltonian. Consider now the following Bryson's formulation:
= h − hmin = 0 S (1) = h& = V sin γ = 0 and V ≠ 0 ⇒ γ (t ) = 0 for t ∈ [t 0 , t1 ] S (2 ) = h&& = V& sin γ + Vγ& cos γ = 0 S
(21)
⇒ γ& (t ) = 0 for t ∈ [t 0 , t1 ]
Thus the constraint is of second order as controls appear in γ& and V& , see equations (10a) and (10b). The Hamiltonian for Bryson's formulation can be defined as: H B = λ B f + µS (2 )
(22)
In contrast to Bryson’s formulation, the Hamiltonian for Jacobson’s formulation is given by H J = λ J f + υS
(23)
Note that λ B ≠ λ J , in general, because of different definitions of the Hamiltonian.
Science and Technology Policy ISTECS Journal Vol. VIII/2006
19
Figure 3: Dircol computational results for the costate variables of Bryson’s formulation.
The direct method approach for optimal control mainly uses Jacobson’s formulation in the derivation of the Karush-Kuhn-Tucker (KKT) necessary conditions. Therefore the costate estimation from the direct method is accurate for the problem having a mixed constraint while it may not work well for the problem having a pure state constraint. Thus, in general, for a pure state constraint situation the direct method will compute λ J , while the indirect method will need λB . The following example gives some insights into the different costate estimation for Bryson’s and Jacobson’s formulation using DIRCOL. 3. NUMERICAL EXAMPLES AND DISCUSSIONS In this example DIRCOL is implemented for the minimum altitude problem where the state constraint is active along the optimal trajectory. Figures 3 and 5 show DIRCOL solutions obtained using the following data:
Science and Technology Policy ISTECS Journal Vol. VIII/2006
20
γ (0) = 0 deg
γ (t f ) = 0 deg
V (0) = 272 m/s
V (t f ) = 306.324004 m/s
x ( 0) = 0 m
x(t f ) = 5813.44774 m
h(0) = 30 m
h(t f ) = 30 m
The following equation shows the differences in the costate equations for Bryson's and Jacobson's formulation. State and costate equations for Bryson's formulation T −D L sin α + cos α − g m m T −D L & cos α − sin α V= ⇒ reduced state equations m m x& = V & h=0
γ& =
λVB ρVS ref [C d cos α + Cl sin α ] + λBx = m Cl ρS ref cos α C d ρS ref sin α T sin α g B − λγ + µV − − 2 + 2 m m 2 2 V m V costate ⇒ =0 equations λVBV 2 S ref ρ h [C d cos α + Cl sin α ] = 2m B V VS λ µ ρ + γ ref h [ − 1 + − C d sin α + Cl cos α ] 2m
λ&γB = λVB g − λhBV − µV& λ&VB
(
λ&Bx λ&hB
)
(
)
Figures 3 and 5 show that the computational results for the costate variables are very different. It is obvious because the constraint which is adjoined to the Hamiltonian differs for both cases. However, the state variables give the same approximate solutions (see Figures 2 and 4). Note that, when solving the TPBVP corresponding to the Jacobson’s formulation, finding υ does not lend itself to a systematic iterative procedure, as pointed out by Maurer and Gillesen [11, Page 111]. On the other hand, µ in the Bryson’s formulation can be found readily using the conditions γ& = 0 and H αB = 0 .
Science and Technology Policy ISTECS Journal Vol. VIII/2006
21
Figure 4: Dircol computational results for the state variables of Jacobson’s formulation.
State and costate equations for Jacobson's formulation T −D L g cos γ =0 γ& = sin α + cos α − mV mV V T −D L V& = cos α − sin α − g sin γ m m x& = V cos γ h& = V sin γ
λ&γJ = λVJ g cos γ − λ&VJ =
λVJ ρVS ref
λγJ V
g sin γ + λ Jx V sin γ − λhJ V cos γ
[C d cos α + Cl sin α ] − λ Jx V cos γ − λhJ V sin γ
m Cl ρS ref cos α C d ρS ref sin α T sin α g − λγJ − − 2 + 2 cos γ 2m 2m V m V
λ&Jx = 0 λ&hJ =
λVJ V 2 S ref ρ h 2m
[C d cos α + Cl sin α ] − υ
λγJ VS ref ρ h [C d sin α − Cl cos α ] − 1 − 2m Science and Technology Policy ISTECS Journal Vol. VIII/2006
22
Figure 5: Dircol computational results for the costate variables of Jacobson’s formulation.
REFERENCES [1] D. Benson. A Gauss Pseudospectral Transcription for Optimal Control. PhD thesis, Massachusetts Institute Technology, February, 2005. [2] A. E. Bryson, W. F. Denham, and S. E. Dreyfus. Optimal programming problems with inequality constraints 1: Necessary conditions for extremal solutions. AIAA Journal, 1(11):2544-2550, 1963. [3] F. Fahroo and I. M. Ross. Costate estimation by a Legendre pseudospectral method. Journal of Guidance, Control, and Dynamics, 24(2):270-277, 2001. [4] W. Grimm and A. Markl. Adjoint estimation from a direct multiple shooting method. Journal of Optimization Theory and Applications, 92(2):263-283, 1997. [5] D. H. Jacobson, M. M. Lele, and J. L. Speyer. New necessary conditions optimality for control problems with state-variable inequality constrains. Journal of Mathematical Analysis and Applications, 35:255-284, 1971. [6] H. Maurer and W. Gillessen. Application of multiple shooting to the numerical solution of optimal control problems with bounded state variables. Computing, 15:105-126, 1975. Science and Technology Policy ISTECS Journal Vol. VIII/2006
23
[7] O. von Stryk. Numerical solution of optimal control problems by direct collocation. In R. Bulirsch, A. Miele, J. Stoer, and K. H. Well, editors, Optimal Control, Calculus of Variations, Optimal Control Theory and Numerical Methods, volume 111 of International Series of Numerical Mathematics, pages 129-143. Birkhäuser Verlag, Basel, 1993. [8] O. von Stryk. User's guide for DIRCOL - a direct collocation method for the numerical solution of optimal control problems. Technische Universität Darmstad, November 1999. [9] O. von Stryk and R. Bulirsch. Direct and indirect methods for trajectory optimization. Annals of Operations Research, 37(1-4):357-373, 1992. [10] O. von Stryk and M. Schlemmer. Optimal control of the industrial robot manutec r3. In R. Bulirsch and D. Kraft, editors, Computational Optimal Control, volume 115 of International Series of Numerical Mathematics, pages 367-382. Birkhäuser, Basel, 1994.
Profile
Subchan received the Bachelor from Institut Teknologi 10 Nopember Surabaya in 1994. He has been worked at the Indonesian Aircraft Industries from 1994-1997. He obtained the Master of Science in Technical Mathematics from Technische Universiteit Delft in 2000 and the Doctorate from College of Defence Technology, Royal Military College of Science, Cranfield University, Defence Academy of the United Kingdom in 2005. Since 2005 he is a research officer in Defence Academy of the United Kingdom. His research interests include dynamic optimization, optimal control, numerical analysis, decision making and applied mathematics.
Science and Technology Policy ISTECS Journal Vol. VIII/2006
24
ISTECS JOURNAL, Vol. VIII (2006) 25 – 33
Fabrication of Bovine Bone Hydroxyapatite: effect of the material shapes and calcination temperature M.K. Herliansyah1, 2, J.A. Toque1, 3, M. Hamdi1, A. Ide-Ektessabi4, M.W. Wildan2 1
Department of Engineering Design and Manufacture, Faculty of Engineering, University of Malaya, Kuala Lumpur,
Malaysia 2Department of Mechanical Engineering, Gadjah Mada University, Yogyakarta, Indonesia 3Department of Mechanical Engineering, University of the Philippines, Diliman, Quezon City, Philippines 4
International Innovation Center, Kyoto University, Sakyo-ku, 606-8501, Kyoto, Japan E-mail:
[email protected]
ABSTRACT The production of hydroxyapatite (HA) from cortical bovine bones was studied in this paper. Bovine hydroxyapatite (BHA) was produced from bovine bone in bulk and powders form by calcination without compaction. The bovine bones were calcined at temperatures ranging from 700-1100ºC. It was discovered that material shape has some influence on the calcination behavior of the bovine bones. XRD results confirmed that HA has been successfully produced from bulk bovine bone but traces of α-TCP and β-TCP were also found in BHA from bovine bone powder. Keywords: bovine bone hydroxyapatite, bulk bovine bone, bovine bone powder, calcination
1. INTRODUCTION In recent years, hydroxyapatite ceramic (HA) implants have attracted attention since it may be possible to use them as an alternative to autogenous free bone grafting. Many basic and clinical studies have confirmed the HA as a potential for biomedical applications, and its clinical applications are gradually expanding [1-7]. Pure HA which is Ca10(PO4)6(OH)2, has chemical compositions closely
resembling that of mineral phase of natural hard tissues of the body such as bone and teeth [6]. It also has high biocompatibility and good bioaffinity, stimulates osteoconduction and is slowly replaces by the host bone after implantation [6]. HA can be used either as a coating material or a bone graft. In coating application, it is normally deposited as a thin-film on to metallic implants using various coating techniques such as PVD, CVD, etc. [5]. As a bone graft, it is being used either as dense or porous grafts. However there are still some problems with this material, such as low porosity and small pore size which is difficult for the bone cells to grow inside, poor interconnection which is important to bone conduction and poor mechanical properties [7]. All these shortcomings limit its applications. Therefore it is a must to find ways how to resolve this issue.
Science and Technology Policy ISTECS Journal Vol. VIII/2006
25
Because of the attractive properties of HA, various techniques have been being developed to produce hydroxyapatite. There are two main ways of producing HA; one is inorganic synthesis such as wet-precipitation, continuous precipitation, sol-gel method, hydrothermal method, thermal deposition, and mechano-chemical method, while another one is from natural source such as corals [8], egg shells [9,10], cuttlefish shells [11, 12], natural gypsum, natural calcite [13], and bovine bone [14]. In this study, bovine bone will be used to prepare HA because it is morphologically and structurally similar to human bone [15]. Furthermore, bovine bone is easy to obtain, lower cost and available in unlimited supply. Two important process parameters i.e. bovine bone raw material shapes and calcination temperature were studied and reported 2. EXPERIMENTAL PROCEDURES 2.1. Sample Preparation Cortical (femur) bovine bones were collected from the local slaughter houses. The procured bone samples were cleaned using boiling method to remove organic substances and collagen. This was done to avoid soot formation in the material during the calcination process. Raw bone was boiled in water for 30 minutes at 99.5oC, and then the water was removed and the bones were washed using fresh water.
This process was repeated three times until it yielded white and clean samples. Before boiling, the macroscopic adhering impurities and substances, which include the ligaments and tissues that stick on the bone were shaved and removed. After boiling, the bone samples were sun-dried for 3 days. The dried cortical bone samples were cut into cubic shape (10 mm x 10 mm x 5 mm) using hacksaw. The bulk and powders bone from cutting chip or sawdust were collected for this experiment. Both of the bovine bone samples in cubic and powder forms were calcined in a box furnace at the following temperatures of 700oC, 800oC, 900oC, 1000oC, and 1100oC; with a heating rate of 5oC/min, in air atmosphere. The respective temperature was maintained for 2 hours to remove the organic matrix. The samples were then cooled to room temperature by slow furnace cooling. 2.2. Characterization The both shapes of HA samples that were obtained from calcined bovine bone at various temperatures were analyzed using XRD with a monochromated CuKα radiation. A scan speed of 7º per minute and a step scan of 0.02º were employed. The chemical elements of the calcined specimens were analyzed using FTIR to confirm the chemical compound after calcination. 3. RESULTS AND DISCUSSION 3.1. General Observation Figure 1 shows the microscopic morphologies of the different shape of bovine bone after calcination at different temperatures. Generally all the samples show a white color and very small size of pores in the surface of bulk samples. Figure 2 shows the percentage of weight reduction of the sample increases with increasing calcination temperature. The highest percentage of weight reduction for bulk samples is 36.1% after calcination at 1100oC for 2 hours and 20% for powder samples after calcination
at 900oC for 2 hours. Science and Technology Policy ISTECS Journal Vol. VIII/2006
26
BULK 700oC
BULK 800oC
POWDER 800oC
POWDER 700oC
BULK 1100oC
BULK 900oC
BULK 1000oC
POWDER 900oC
POWDER 1000oC
o POWDER 1100 C
Figure 1. Appearances of bovine bone sintered at different shapes (bulk and powder) and various heating temperature
Bulk
40.00 % Weight Reduction
Powder
35.00 30.00 25.00 20.00 15.00 600
700
800
900
1000
1100
1200
Calcine Temperature (Degree Celcius)
Figure 2. Effect of calcination temperature to samples weight reduction
The average of weight reduction (34.53%) for bulk samples is higher than of powder samples (19.2%). This different weight reduction is due to the organic content in the bulk samples higher than that of in powder samples. This result corresponds to XRD analysis that shown the decomposition phenomena from HA to β-TCP and α-TCP. 3.2. XRD Phase Analysis The results of the X-ray diffraction experiments that are indicative for the chemical composition (presence of crystalline phases) are shown in Figure 3 (a). for BHA from bulk bovine bone and Figure 3 (b). for BHA from bovine bone powder. The XRD pattern of the specimens showed peaks characteristic of that of HA (JCPDS 9-432). All of the peaks in Figure 3 (a). belong to HA, and they all exhibit a high Science and Technology Policy ISTECS Journal Vol. VIII/2006
27
crystallinity as indicated by the narrow diffraction peaks. The increasing of cristallinity together with increasing of calcinations temperature was also shown in Figure 3 (a).
a
b
* = α-TCP
x = β-TCP
Figure 3. XRD Result of BHA from calcined (a) granule and (b) powder bovine bone
However the result of the characterization of the BHA from bovine bone powder has shown some intriguing result. The XRD patterns, as shown in Figure 3 (b). of the bovine bone powder calcined at different temperatures ranging from 700-1100ºC show small amounts of impurities traces of α-TCP and β-TCP besides the major phase (peaks marked with asterisks and x in Fig. 3(b)) even when it was calcined at relatively low temperatures. This means that the decomposition process of HA to other phases occured at a lower temperature of what is expected. Previous studies showed that such phase transformation will start at temperature above 1000oC [14, 16-18]. In this research, these phase were detected to have been present with temperature as low as 700oC during calcinations of bovine bone in the form of powder. One of the reasons would be material shape. The previous result (Figure 3 (a).) used same material (bovine bone) but different shape. The result in Figure 3 (b). used bovine bone powders which are not compacted. Therefore these powders have significant amount of porosity in it to allow the transfer of heat more effectively. This owes to the fact that powders have bigger surface area to volume ratio. This may not seem significant on a large scale but since powders are on the scale of microns to sub-microns, this has tremendous effect. The ratio could go up to an order of as high as 106-109. Having a large ratio can mean that the effect of heating can occur at much lower expected temperature. It should be noted that the bovine powders were needed to be calcined to produce the HA as compared to the commercial HA where the heating process was performed to sinter the compacted powders not to produce it. In addition, since the specimens were in powder form, the organic compounds found in the bulk bone can be easily removed and therefore will require lower temperature and lesser time to initiate the growth of HA phase. The XRD result of BHA from bulk bovine bone in Figure 3 (a). supports the argument presented. It shows the XRD pattern of bulk bovine bone calcined at the same conditions from Science and Technology Policy ISTECS Journal Vol. VIII/2006
28
a previous study [14]. It can be seen that the transformation of HA to β-TCP for the bulk samples in Figure 3 (a). did not start to occur at a calcination temperature of 1100ºC which is in agreement to other previous studies. The experiment of Figure 3 (a). and Figure 3 (b). only differs in the material shape. Both of the samples were calcined after defatting and got from the same bovine bone. Moreover, Figure 3 (a). shows the XRD pattern of the raw bovine bone. As expected, it resembles a pattern of that of an amorphous phase. It can be seen that the peaks which manifested after subsequent heating was not present from the starting material and that the calcium phosphate phases identified after heating resulted from the calcination process. In addition, the presence of other calcium phosphate phases may not be all bad. Some recent studies revealed that biphasic calcium phosphate produced good results in biological studies [19-25], while other studies showed that the presence of these other phases contribute to the rapid dissolution when implanted in to the body [25]. These phases coexisting with HA might provide better anchorage for bony tissues since they dissolve faster in physiological fluid. The next step of this study is to evaluate the mechanical properties and biological properties of the BHA 3.3. FTIR Analysis The chemical compounds present in the BHA samples were determined by FTIR spectroscopy. The FTIR spectrum of BHA obtained from KBr pellets is shown in Figure 4 for bulk bovine bone (a) and bovine bone powder (b) calcined at temperatures of 700oC, 800oC, 900oC, 1000oC, and 1100oC
respectively. Both of the IR spectrum exhibit only the characteristic absorption peaks of HA. A large number of bands in the spectra matches the bands in the HAp reference spectrum and are in agreement with literature data on HA [26, 27].
a
b
Figure 4. FTIR result of BHA from calcined bovine bone (a) bulk and (b) powder Science and Technology Policy ISTECS Journal Vol. VIII/2006
29
60
OH CO
50
PO
30
CO
25 Intensity
40 Intensity
OH
35
PO
30
20 15
20
10
10
5
0 600
700
800
900
1000
Temperature (degree)
1100
1200
0 600
700
800
900
1000
1100
1200
Temperature (degree Celcius)
a b Figure 5. Effect of calcination temperature and intensity of BHA compound from bovine bone (a) bulk and (b) powder
Figure 4 and Figure 5 also shown that OH- compound and PO43- are decrease related with the increasing of calcinations temperature. It reveals a strong absorption of CO32- group at 1400 -1570 which is lattice site of PO43, but CO32- also decrease related with the increasing of calcinations temperature. 4. CONCLUSIONS Bovine hydroxyapatite was successfully produced by calcinating bovine bone in bulk and powders form. The results showed that BHA can be harnessed at a much lower temperature when calcined as a powder without compaction as compared to calcination of the bulk bones. This finding may be quite promising since this will lessen the time of BHA preparation and will also enable it to be produced at lower temperatures. It was also found out that calcinating in powder form lowers the temperature by which α-TCP and β-TCP start to decompose due to the higher surface to volume ratio of powders as compared to bulk form. The next step for this research is to study the Ca/P ratio and mechanical properties as well as the biological properties of the BHA produced. Low temperature calcination may also be studied. ACKNOWLEDGMENT This work was supported by a grant from the JICA-AUN/SEED Net collaborative research program, Department of Engineering Design and Manufacture, University of Malaya, Malaysia, and Department of Mechanical and Industrial Engineering Gadjah Mada University, Indonesia. REFERENCES [1] Greenspan, D. C., (1999), Bioactive ceramic implant materials, Current Opinion in Solid State and Materials Science 4 pp 389–393. [2] FinnRA, Bell WH, Brammer JA. Interpositional grafting with autogenous bone and corlline hydroxyapatite. J Maxillofac Surg 1980;8:217-227 [3] Shikinami, Y. and Okuno, M, (1999), Bioresorbable Devices Made of Forged Composites of
Hydroxyapatite (HA) Particles and Poly-L-Lactide (PLLA): Part I. Basic Characteristics, Science and Technology Policy ISTECS Journal Vol. VIII/2006
30
Biomaterials 20 pp 859-877 [4] Tadic, D., Epple, M., (2003), Mechanically Stable Implants Of Synthetic Bone Mineral By Cold Isostatic Pressing, Biomaterials 24 pp 4565–4571 [5] Toque J, Hamdi M, Ide-Ektessabi A (2006) A review on hydroxyapatite coating using magnetron sputtering, ICMM Proc. vol. 1, 1st International Conference & 7th AUN/SEED-Net Fieldwise Seminar on Manufacturing and Material Processing, Kuala Lumpur, Malaysia, 2006, pp 603-608. [7] Hedia, H.S. and Nemat-Alla Mahmoud, (2004) Design optimization of functionally graded dental implant, Bio-Medical Materials and Engineering 14 (2004) 133–143 133, IOS Press [8] Herliansyah M. K., Hamdi M, Ide-Ektessabi A, Wildan M. W. (2006). Fabrication of hydroxyapatite bone graft for implant application: a literature study, ICMM Proc. vol. 1, 1st International Conference & 7th AUN/SEED-Net Fieldwise Seminar on Manufacturing and Material Processing, Kuala Lumpur, Malaysia, 2006, pp 559-564. [9] Murugan, R., and Rao, K. P., (2002), Controlled Release of Antibiotic from Surface Modified Coralline Hydroxyapatite, Trends Biomater. Artif. Organs. Vol 16 (1) pp 43-45 [10] Sasikumar, S. and Vijayaraghavan, R., (2006), Low Temperature Synthesis of Nanocrystalline Hydroxyapatite from Egg Shells by Combustion Method Trends, Biomater. Artif. Organs, Vol 19(2),pp 70-73. [11] Prabakaran, K., Balamurugan, A., And S Rajeswari, S., (2005), Development Of Calcium Phosphate Based Apatite From Hen’s Eggshell, Bull. Mater. Sci., Vol. 28, No. 2, pp. 115–119, Indian Academy of Sciences. [12] Rocha, J.H.G., Lemos, A.F., Agathopoulos, S., Valério, P., Kannan, S., Oktar, F.N., and Ferreira, J.M.F., (2005), Scaffolds for bone restoration from cuttlefish, Bone 37 pp 850–857 [13] Rocha, J.H.G., Lemos, A.F., Kannan, S., Agathopoulos, and Ferreira, J.M.F., (2005), Hydroxyapatite Scaffolds Hydrothermally Grown from Aragonitic Cuttlefish Bones, J. Mater. Chem., 15, pp 5007–5011 [14] Nasution, D. A., (2006), Fabrication and Analysis of Physical and Mechanical Properties of Hydroxyapatite from Natural Calcite, Master Thesis, Department of Mechanical Engineering, Faculty of Engineering, Gadjah Mada University, Yogyakarta, Indonesia [15] Ooi, C. Y., M. Hamdi, M., and Ramesh, S., (2006), Properties of hydroxyapatite produced by annealing of bovine bone, Ceramics International, doi:10.1016/j.ceramint. 2006.04.001. [16] Goren, S., Gokbayrak, H. and Altinas, S. (2004). Production of Hydroxyapatite from Animal Bone. Key Engineering Materials, 264-268, 1949-1952 [17] Muralithran G, Ramesh S (2000). The effects of sintering temperature on the properties of hydroxyapatite. Ceram Int 26: 221-230. [18] Tampieri A, Celotti G, Szontagh F, Landi E (1997). Sintering and characterization of HA and TCP bioceramics with control of their strength and phase purity. J Mat Sci: Mat Med 8: 29-37. [19] Rao W and Boehm R (1974). A study of sintered apatites. J Dent Res 53: 1351-1354. [21] Kalita S, Bhardwaj A, Bhatt A (2006). Nanocrystalline calcium phosphate ceramics in biomedical engineering. Mat Eng C (In Press, Corrected Proof). Science and Technology Policy ISTECS Journal Vol. VIII/2006
31
[22] Trojani C, Boukhechba F et al (2006). Ectopic bone formation using an injectable biphasic calcium phosphate/Si-HPMC hydrogel composite loaded with undifferentiated bone marrow stromal cells. Biomaterials 27: 3256-3264. [23] Nihouannen D, Guehennec L et al (2006). Micro-architecture of calcium phosphate granules and fibrin glue composites for bone tissue engineering. Biomaterials 27: 2716-2722. [24] Santos E, Farina M, Soares G (2006). Specific proliferation rates of human osteoblasts on calcium phosphate surfaces with variable concentrations of α-TCP. Mat Sci & Eng C (In Press, Corrected Proof). [25] Goyenvalle E, Aguado E et al (2006) .Osteointegration of femoral stem prostheses with a bilayered calcium phosphate coating. Biomaterials 27:1119-1128. [27] Curran J, Gallagher J, Hunt J (2005).The inflammatory potential of biphasic calcium phosphate granules in osteoblast/macrophage co-culture. Biomaterials 26: 5313-5320. [28] Yang Y, Kim K-H, Ong J (2005). A review on calcium phosphate coatings produced using a sputtering process—an alternative to plasma spraying. Biomaterials 26: 327-337. [19] Murugan, R., Sampath Kumar, T. S. and Rao, P. K. (2002). Fluorinated Bovine Hydroxyapatite: Preparation and Characterization. Materials Letters, 57, 429- 433. [20] Joschek, S., Nies, B., Krotz, R. and Gofpferich, A. (2000). Chemical and Physicochemical Characterization of Porous Hydroxyapatite Ceramics Made of Natural Bone. Biomaterials, 21, 1645-1658.
Science and Technology Policy ISTECS Journal Vol. VIII/2006
32
Muhammad Kusumawan Herliansyah was born in Kulon Progo (Indonesia) in 1971. In 1996 he obtained his undergraduate in Mechanical Engineering from the Gadjah Mada University (Yogyakarta, Indonesia). In 2002 he obtained his master degree in Industrial Engineering from Institute of Teknologi Bandung (Bandung, Indonesia). Since December 2005, he has been a PhD student at Department of Engineering Design and Manufacture, University of Malaya under the supervision of Dr. Mohd. Hamdi Abdul Shukor. The objective of his research has been the development and fabrication of hydroxyapatite bone graft for implant application. This research work has also involved other institutional partners, namely the International Innovation Center of Kyoto University (local supervision: Prof. Dr. Ari Ide-Ektessabi) and Department of Mechanical Engineering of Gadjah Mada University (local supervision: Dr. Muhammad Waziz Wildan). Jay Arre O Toque was born in Cavinti, Laguna (Philippines) in 1981. He graduated with a degree in Mechanical Engineering from the University of the Philippines (Diliman, Quezon City, Philippines). Since June 2005, he has been a Master student at Department of Engineering Design and Manufacture, University of Malaya under the supervision of Dr. Mohd. Hamdi Abdul Shukor. His research topic is fabrication of hydroxyapatite coating for biomedical implant application using magnetron sputtering. Mohd. Hamdi bin Abdul Shukor was born in Muar (Malaysia). In 1994 he obtained his B.Eng (Hons) in Mechanical Engineering from Imperial College of Science, Technology and Medicine, University of London, London, United Kingdom. In 1996 he obtained his M. Sc. Degree in Advanced Manufacturing Technology and Systems Management from University of Manchester Institute of Science and Technology (UMIST), Manchester, United Kingdom. He is received the Dr. Eng. in Preparation and characterization of Calcium Phosphate coating for Biomedical applications from Kyoto University (Kyoto, Japan). He currently has a position as a senior lecturer at the Department of Engineering Design and Manufacture, University of Malaya (Kuala Lumpur, Malaysia). He is currently author of more than 9 papers in international refereed journals. Ari Ide-Ektessabi currently has a position as Professor in the International Innovation Center of Kyoto University. His areas of research interest are now in bio Microsystems Electronics, Applications of Synchrotron Radiations to Micro Biology, and Ion Beam Processing and Analysis.
Muhammad Waziz Wildan was born in 1968 in Boyolali, Indonesia. He was graduated in Mechanical Engineering from the Gadjah Mada University (Yogyakarta, Indonesia). In 1997 he obtained his Master Degree in Energy System and Environments – Mechanical Engineering from University of Strathclyde (Glasgow, UK). He is received the Ph.D. in Engineering Materials (Advanced Ceramics) – Mechanical Engineering from the University of Strathclyde (Glasgow, UK) in 2000. He currently has a position as a senior lecturer at the Department of Mechanical Engineering, Gadjah Mada University. As a lecturer, he is responsible for Advanced Ceramics Materials. Science and Technology Policy ISTECS Journal Vol. VIII/2006
33
ISTECS JOURNAL, Vol. VIII (2006) 34 – 43
PID Parameters Optimization by Using Genetic Algorithm A Study on Time-delay Systems Andri Mirzal, Masashi Furukawa Graduate School of Information Science and Technology Hokkaido University Sapporo, Japan Email: andri,
[email protected]
ABSTRACT Time delays are components that make time-lag in systems response. They arise in physical, chemical, biological and economic systems, as well as in the process of measurement and computation. In this work, we implement Genetic Algorithm (GA) in determining PID controller parameters to compensate the delay in First Order Lag plus Time Delay (FOLPD) and compare the results with Iterative Method and Ziegler-Nichols rule results. Keywords: PID controller, time delay, Genetic Algorithm, Iterative Method, Ziegler-Nichols rule
1. INTRODUCTION In the previous work [1], authors have implemented and compared two tuning methods, Iterative Method and Ziegler-Nichols rule, to compensate the effect of delay in stability of systems and showed that Iterative Method has superior performance in analyzed cases, FOLPD (First Order Lag plus Time Delay). But there are some cases where we can’t use these two tuning methods, i.e. the dynamic plants which its parameters are constantly changing. In this sort of systems, we have to do retuning in real time, which can’t be applied by the tuning methods because we have to take the systems offline first in order to set its parameters [3].
1 Ti s R(s) +
+
E(s)
X
Kc
-
+ U(s)
X
G p ( s )e
− sτ p
Y(s)
+
Td s
Figure 1. The PID controller general structure where plant has delay component Science and Technology Policy ISTECS Journal Vol. VIII/2006
34
In this work, we extend our previous work [1] by implementing Genetic Algorithm (GA) in determining PID Controller parameters to compensate the delay in order one (FOLPD) and compare the results with Iterative Method and Ziegler-Nichols rule results. 2. THE OBJECTIVE FUNCTIONS FITNESS VALUES The most crucial step in applying GA is to choose the objective functions that are used to evaluate fitness of each chromosome. Some works [3] [4] use performance indices as the objective functions. In [3] author uses Mean of the Squared Error (MSE), Integral of Time multiplied by Absolute Error (ITAE), Integral of Absolute Magnitude of the Error (IAE), and Integral of the Squared Error (ISE), while in [4] authors use ISE, IAE, and ITAE. Here we use all four performance indices stated above and Integral of Time multiplied by the Squared Error (ITSE) to minimize the error signal E(s) and compare them to find the most suitable one. The performance indices are defined as follow [2]: τ
τ
τ
τ
τ
1 2 MSE = ∫ (e(t ) ) dt , ITAE = ∫ t e(t ) dt , IAE = ∫ e(t ) dt ISE = ∫ e(t ) 2 dt , and ITSE = ∫ te(t ) 2 dt t0 0 0 0 0
(1)
Where e(t) is the error signal in time domain. The PID controller is used to minimize the error signals, or we can define more rigorously, in the term of error criteria: to minimize the value of performance indices mentioned above. And because the smaller the value of performance indices of the corresponding chromosomes the fitter the chromosomes will be, and vice versa, we define the fitness of the chromosomes as:
fitness value =
1 performance index
(2)
3. DELAY COMPONENT Delay in control systems can be defined as time-interval between an event that start in one point with its output in another point within systems [5]. Delay is also recognized as transport lag, deadtime, and time lag. Because delay always reduces stability of minimum phase systems (systems which don’t have poles and zeros in the right half of s-plane), it is important to analysis stability of systems with time delay. In order to do tuning processes using GA, we approximate the delay with Direct Frequency Response (DFR) series. because it has been shown in [1] that this series has the smallest average error among the others seven series. Furthermore, for simulation purpose, we use second order DFR series to circumvent unnecessary complexity, because as the order of the series are getting higher, not only the calculation becomes difficult but also it introduces new poles and zeros which make the system much more elusive.
DFR series : Science and Technology Policy ISTECS Journal Vol. VIII/2006
e − sτ ≈
1 − 0.49sτ + 0.0954s 2τ 2 1 + 0.49sτ + 0.0954s 2τ 2
(3) 35
4. GENETIC ALGORITHM PID controller parameters will be optimized by applying GA. Here we use Matlab Genetic Algorithm Toolbox [6] to simulate it. The first and the most crucial step is to encoding the problem into suitable GA chromosomes and then construct the population. Some works recommend 20 to 100 chromosomes in one population. The more the chromosomes number, the better the chance to get the optimal results. However, because we have to consider the execution time, we use 80 or 100 chromosomes in each generation. Encoding is done in real number rather than binary encoding because the latter discards the parameters value if it exceeds its precision capability. Each chromosome comprises of three parameters, Kd, Kp, Ki, with value bounds varied depend on the delay and objective functions used. After many
experiments, we find that the value bounds should be set according to the Iterative Method and Ziegler-Nichols rule value range to ensure the convergence (there are many cases which the convergence can’t be reached if we set the parameter value bounds arbitrarily, even though the optimal results included in those bounds range). The population in each generation is represented by 80 x 4 or 100 x 4 matrix, depends on the chromosomes number in population, which each row is one chromosome that comprise Kd, Kp, Ki values and the last column added to accommodate fitness values (F) of corresponding chromosomes. K d1 K d2 . . K dn
K p1
K i1
K p2
K i1
.
.
. K pn
. K in
F1 F2 . . Fn
chromosome1 chromosome 2 .
(4)
. chromosome n
We use maximum generation termination (maxGenTerm.m) to terminate the program rather than considering the best chromosome fitness values changing rate because we want to control the execution time. However, the best chromosome fitness values changing rate is also being considered by running the programs until the best fitness value stop increasing, then we set that point as the maximum generation. After several experiments it’s shown that there is no visible improvement after 300th generation, so we set 300 as the maximum generation. Matlab GA Toolbox [6] provides three selection techniques, Tournament Selection, Roulette Wheel Selection and Normalize Geometric Selection. Tournament Selection requires more execution time while Roulette Wheel Selection allows the weaker chromosomes to be selected many times, so we choose Normalized Geometric Selection to choose the parents. After parents being selected, the crossover operation will be done. We use arithmetic crossover (arithXover.m) function because it is specifically being used for floating point numbers and provides more than one crossover points. And we set four crossover points because our chromosome comprise of three alleles, one point crossover can not accommodate three alleles in one operation. Mutation is done by setting mutation probability around 0.1 percent. In general mutation operations should not be done too often because the searching process will change into random search as Science and Technology Policy ISTECS Journal Vol. VIII/2006
36
the mutation probability getting higher. 5. APPLYING THE GENETIC ALGORITHM There are several variables used as the standard to measure systems performance. In general, unit step input is used to test the systems, and the output signals is characterized by some standard performance measures: settling time (ST), percent overshoot (PO), error signal, rise time (RT), peak time (PT), and stability margin SM). All these measures are defined in time domain response. 5
60 Ziegler Nichols Iterative Method MSE IAE ISE ITAE ITSE
(a) 50
(b)
4.5
4
3.5 40
t o o h sr e v O t n e cr e P
n oi r eti r c % 5 e m ti g nli tt e S
30
20
3
2.5
2
1.5
Ziegler Nichols Iterative Method MSE IAE ISE ITAE ITSE
1 10
0.5
0
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
0
1
0
0.1
0.2
0.3
0.4
1
0.7
0.8
0.9
1
0.6
0.7
0.8
0.9
1
10 Ziegler Nichols Iterative Method MSE IAE ISE ITAE ITSE
(c)
Ziegler Nichols Iterative Method MSE IAE ISE ITAE ITSE
(d) ) el a c s ci m hti r a g o(l d n o c e s ni e m ti k a e P
0
10
-1
10
1
10
0
10
-1
10
-2
-2
10
0.6
2
10
) el a c s ci m hti r a g o(l d n o c e s ni e m ti e si R
0.5
Delay in second
Delay in second
10 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
0
0.1
0.2
0.3
0.4
0.5
Delay in second
Delay in second 3
10
) el a c s ci m hti r a g o(l ni g r a m yt ili b at S
Ziegler Nichols Iterative Method MSE IAE ISE ITAE ITSE
(e) 2
10
1
10
0
10
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Delay in second
Figure 2. Comparison of Standard performance Measures, (a) PO, (b) ST, (c) RT, (d) PT, (e) SM.
Science and Technology Policy ISTECS Journal Vol. VIII/2006
37
Figure 2 describes standard performance measures of a typical system driven by unit step input. Percent overshoot is defined as the point where the system response reaches the peak, in this case 53%. There are several criteria for settling time, for example 1% criterion, 2% criterion, and 5% criterion. Here we use 5% criterion settling time. And for the rise time, actually in general, is measured as the time needed by systems to reach from 0 to 100% of final value or from 10% to 90% of final value. But, for measurement simplicity, we use 0 - 95% criterion. Peak time is the point where the maximum value reached. And error signal is the difference between the input signal magnitude and system response final magnitude. In this work, we use G(s) = 1/s+1, delay is in the range of 0.01 to 1 second. And because the systems are compensated by PID controller, the error signals are always zero. In addition to the five system standard performance measures described above, in Iterative Method and Ziegler-Nichols rule, we calculate the system performance indices described by equation (1) also. This is done because we want to compare it with the result of GA, which it’s being optimized in the term of performance indices. Ideally, we can expect corresponding GA’s performance indices should be always better than two tuning rules. To calculate performance indices, we approximate the integral in equation (1) with addition (sigma) and 0.01 second sampling time and set the sigma upper limit with 15 second for all analyzed cases, no matter how quick it reaches convergence values. 6. SIMULATION RESULTS AND ANALYSIS Table 1. Average values of standard performance measures. Parameter
Ziegler-Nichols
Iterative
Optimized
Optimized
Optimized
Optimized
Optimized
Method
by MSE
by IAE
by ISE
by ITAE
by ITSE
%OV
38%
15%
10%
6%
11%
10%
8%
ST5%
1.53
1.54
1.47
1.01
1.45
0.745
1.37
RT
0.444
0.912
0.453
0.495
0.455
0.588
0.474
PT
0.613
3.43
0.576
0.622
0.576
0.836
0.597
SM
36.25
37.86
33.68
33.8
24.4
36.8
32
6.1. Standard performance measures.
Percent Overshoot While figure 2(a) summarizes the values change of percent overshoots with respect to the time delay, table 1 gives its average values. Ziegler-Nichols rule gives the biggest value for all time delay, consequently its average value is the biggest also, 38%. Here the difference between two tuning methods and GA methods can be seen: while tuning methods have almost the same pattern, its value decreasing as the time delay growing bigger, except for small value of delay where Ziegler-Nichols rule gives Science and Technology Policy ISTECS Journal Vol. VIII/2006
38
increasing values, the GA methods give almost constant value, around 10% as long the time delay is not too small, except for case optimized by ITAE, where the percent overshoot value fluctuates and decreasing as time delay increasing. Table 1 shows that GA produces much better percent overshoot than two others tuning methods, especially if optimized by IAE criterion. So it can be concluded that GA can be used to optimized percent overshoot. Settling Time The value of settling time (5% criterion) all over the time delay is summarized by figure 2(b), where it can be seen that almost all methods, except for GA optimized by IAE and ITAE, fall under almost the same straight line with positive slope. It means that as the delay increasing, the settling time will increase linearly. In table 1 we can see that the settling time average value is not so different among all methods, in the range 1.37 second to 1.54 second, except for GA optimized by ITAE 0.745 second and IAE 1.01 second. However, because the others GA methods results are in agree with two tuning methods, we can’t really differentiate these results. So it can be said that the settling time is not optimized by GA methods. Rise Time The third variable is rise time, which plotted in figure 2(c) using logarithmic scale in y axis. We can see strong pattern which all the results, except for Iterative Method, have almost the same value all over the time delay. And it should not be surprising if we get almost the same average value for all result, around 0.44 second to 0.59 second, except for Iterative Method, 0.912 second. Another interesting thing is, in general, almost in all range of time delay the curves keep their ranking unchanged, with order from the biggest value: Iterative Method, ITAE, IAE, ITSE, ISE, MSE, and Ziegler-Nichols rule. From the average value (table 1), the best result is given by Ziegler-Nichols rule, 0.444 second and the worst one is Iterative Method 0.912 second, and all GA methods produce almost the same average value. But because GA results are not so different compared to Ziegler-Nichols rule, it can’t be concluded that GA can optimize the rise time. Peak Time Almost the same pattern, as in rise time plots, is shown on the peak time plots on figure 2(d), except Iterative Method, in large delay range, tends to diverge, where all methods show almost the same values for all over time delay. But we must pay attention on GA optimized by ITAE because there is a range which its value is bigger than the others. Except for Iterative Method which is 3.43 second, all others methods produce peak time around 0.57 second to 0.83 second. The best values are given by GA optimized by MSE and ISE, 0.576 second. Like rise time, in general, almost in all range of time delay the curves keep their ranking unchanged, with order from the biggest value: Iterative Method, ITAE, IAE, Ziegler-Nichols, ITSE, MSE, and ISE. And because the GA methods produce peak time plots which only better than Iterative Method, not Ziegler-Nichols, the peak time can’t be optimized by GA methods. Science and Technology Policy ISTECS Journal Vol. VIII/2006
39
Stability Margin The last standard performance measure is stability margin (figure 2(e)). Stability margin is the maximum gain that can be set before system response goes into sinusoidal cycle. In the simulations, this is done by simply increasing value of Kc until sinusoidal cycle happens, and the stability margin of corresponding system is Kc at sinusoidal cycle. This is the first result that shows consistency all over time delay range, which all curves fall under almost the same line. So this is the strongest patterns, and because the stronger the pattern, the less the ability of GA methods to optimize corresponding performance measures, we can’t use GA methods to optimize the stability margin. Conversely, the GA methods are likely to produce less stable systems. Like rise time and peak time, in general, almost in all range of time delay the curves keep their ranking unchanged, with order from the biggest value: Iterative Method, ITAE, Ziegler-Nichols, IAE, ITSE, ISE, and MSE. But unlike settling time, rise time, and peak time, stability margin values reduce as the time delay increases. 6.2. Performance Indices We differentiate the term standard performance measures and performance indices here, where standard performance measures have already been discussed above, performance indices are: MSE, IAE, ISE, ITAE, and ITSE. Table 2 below summarizes performance indices from simulation results. As expected, the average values of GA performance indices are always smaller than its corresponding Ziegler-Nichols and Iterative Method. Moreover Ziegler-Nichols rule produces smaller average performance indices values than Iterative Method does for all time delay values range. Table 2. Performance indices of tuning methods and GA methods. Ziegler-Nichols
Iterative Method
Optimized
Optimized
Optimized
Optimized
Optimized
Delay MSE
IAE
ISE
ITAE
ITSE
MSE
IAE
ISE
ITAE
ITSE
by MSE
by IAE
by ISE
by ITAE
by ITSE
0.01
0.000973
2.992285
1.461115
0.073894
0.020615
0.000927
2.730597
1.391095
0.061045
0.014833
0.00081
1.623776
0.840777
0.011395
0.005388
0.025
0.002626
6.927413
3.942144
0.432876
0.124957
0.002313
6.300772
3.471926
0.391845
0.089213
0.001884
3.701027
2.827366
0.218351
0.038659
0.05
0.004988
12.98842
7.487622
1.580628
0.452425
0.004376
11.6831
6.567769
1.452929
0.318439
0.003708
7.09663
5.57464
0.301924
0.155625
0.075
0.007197
18.53609
10.80309
3.249557
0.929874
0.006334
16.58511
9.507315
2.979748
0.642704
0.005508
10.75044
8.266998
0.754851
0.350289
0.1
0.009284
23.69959
13.93484
5.304512
1.517501
0.008226
20.99624
12.34746
4.801577
1.035597
0.0073
13.86607
10.95693
1.526731
0.622276
0.25
0.020422
48.31448
30.65349
20.82159
6.345217
0.019092
40.69953
28.65714
15.58065
4.507248
0.017983
34.18455
26.99306
7.69427
3.854827
0.5
0.037119
76.26794
55.71523
45.56416
17.73724
0.037603
71.60299
56.44149
37.68745
16.79849
0.035612
67.45696
53.45713
30.00848
15.15463
0.75
0.05416
108.6372
81.29453
91.98153
36.10692
0.057392
123.8553
86.14606
146.8766
42.31687
0.053071
100.0041
79.65463
63.58332
33.54903
1
0.072117
146.0676
108.247
174.114
65.18471
0.077665
176.087
116.5754
298.2008
81.12556
0.070379
134.0671
105.6145
111.3298
58.69364
Av.
0.02321
49.38123
34.83767
38.12475
14.26883
0.02377
52.28229
35.6784
56.44808
16.31655
0.021806
41.41674
32.68734
23.93657
12.4916
Although it can be seen that MSE has the smallest average values for all three methods, and IAE has the largest average values for Ziegler-Nichols and GA methods, and only second in Iterative Science and Technology Policy ISTECS Journal Vol. VIII/2006
40
Method, it doesn’t imply that one must use MSE and must avoid using IAE as a objective function in GA because these performance indices, as shown by equations (1), have different definitions and cannot be compared. Moreover, as shown in table 1, GA method optimized by MSE doesn’t produce the best results for all analyzed standard performance measures, only for peak time. To get more insight about the comparison among these methods, we plot values change of each performance indices with respect to time delay below. 180
0.08
GA Optimized by IAE IAE of Iterative Method IAE of Ziegler-Nichols rule
160
0.07
140 0.06 120 0.05
e ul a v E S M
e ul a v E AI
(a). The comparison of MSE values. 0.04
100
(b). The comparison of IAE values.
80
0.03 60
0.02 40
0.01
20
GA Optimized by MSE MSE of Iterative Method MSE of Ziegler-Nichols 0
0 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Delay in second
Delay in second
(a). The comparison of MSE values.
(b). The comparison of IAE values.
120
300 GA Optimized by ITAE ITAE of Iterative Method ITAE of Ziegler-Nichols
e ul a v E SI
100
250
80
200
e ul a v E A TI
60
150
40
100
20
50 GA Optimized by ISE ISE of Iterative Method ISE of Ziegler-Nichols
0
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
0
1
0
0.1
0.2
0.3
0.4
Delay in second
0.5
0.6
0.7
0.8
0.9
1
Delay in second
(c). The comparison of ISE values.
(d). The comparison of ITAE values.
0.08
0.07
0.06
0.05
e ul a v E S TI
0.04
0.03
0.02
0.01 GA Optimized by ITSE ITSE of Iterative Method ITSE of Ziegler-Nichols 0
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Delay in second
(e). The comparison of ITSE values Figure 3. Performance Indices Science and Technology Policy ISTECS Journal Vol. VIII/2006
41
All five figures confidently show that GA method gives the smallest values of all performance indices analyzed for all range of time delay. So not only in average values, but also for all measured values does GA method produce the smallest corresponding performance indices. However, the differences between GA method and two tuning methods results, except for ITAE objective function where the differences increase as the time delay increases, are not impressive enough to come into conclusion that GA method is much better than two others methods in minimizing error criteria. Besides, we must consider the convergence problem arises in applying GA, which in this work the experiments don’t always come into desired solutions. Even though we set the value bounds based on the previous results from two tuning methods, it only improves the probability that the simulations come into convergence results (from 45 cases, there two times failed to reach convergence results).
7. CONCLUDING REMARKS Genetic Algorithm applied in PID controller improves FOLPD transient response compared to two tuning methods. This is shown by average percent overshoot reduction, more than 70% and 30% with respect to the Ziegler-Nichols rule and Iterative Method, while keep the rise time and peak time almost unchanged and improves the settling time. However, there is payoff in the stability margin which reduces slightly compared to two tuning methods. The average values of GA performance indices, as expected, are always smaller than its corresponding Ziegler-Nichols and Iterative Method. Moreover Ziegler-Nichols rule produces smaller average performance indices values than Iterative Method does for all time delay values range. However, the differences between GA method and two tuning methods results, except for ITAE objective function where the differences increase as the time delay increases, are not impressive enough to come into conclusion that GA method is much better than two others methods in minimizing error criteria. There are convergence problems that arise in applying GA, which in this work the experiments don’t always come into desired solutions. Even though we set the value bounds based on the previous results from two tuning methods, it only improves the probability that the simulations come into convergence results. Moreover, the value bounds setting based on tuning methods results discards the possibility to find optimum results from others value ranges. REFERENCES
[1] Andri Mirzal, Shinichiro Yoshii, Masashi Furukawa, Approximation and Compensation of Delay in Analog Control Systems, 精密工学会北海度支部学術講演会, Sapporo, Japan, 2006. [2] Richard C. Dorf, Robert H. Bishop, Modern Control Systems 10th Edition, Pearson Prentice Hall, 2005. [3] Ian Griffin, On-line PID Controller Tuning using Genetic Algorithms, Dublin City University, 2003. [4] T. O’Mahony & CJ. Downing (Cork Institute of Technology, Ireland), Klaudiusz Fatla (Wroclaw University of Technology, Poland), Genetic Algorithms for PID Parameter Optimization, Science and Technology Policy ISTECS Journal Vol. VIII/2006
42
Minimizing Error Criteria. [5] O’ Dwyer, A., The Estimation and Compensation of Processes with Time Delays, Ph.D. Thesis, School of Electronic Engineering, Dublin City University, 1996. [6] C. R. Houck, J. Joines. and M.Kay. A Genetic Algorithm for Function Optimization: A Matlab Implementation. ACM Transactions on Mathematical Software, 1996.
Profile Andri Mirzal received his Bachelor degree (Honor) in Electrical Engineering majoring in Control Engineering from Bandung Institute of Technology, Indonesia, 2003. He is currently a student in the Graduate School of Information Science and Technology. His current interests are in optimization methods, control systems, and network theories. He is now carrying out his research program under supervision of Professor Masashi Furukawa, Autonomous Systems Engineering Laboratory, Graduate School of Information Science and Technology, Hokkaido University. Masashi Furukawa received the B.Sc., M.Sc., and Ph.D. degrees from the University of Hokkaido, Sapporo, Japan. He is currently a Professor in the Graduate School of Information Science and Technology, Hokkaido University, Sapporo, Japan. During 1970s, he was engaged in the development of 3-D CAD/CAM systems. He was a Research Associate participating in the Cornell Injection Molding Project in 1976 – 1977 and in the Computational Geometry Project of the University of East Anglia in 1981 – 1982. His current interests are machine learning and GAs.
Science and Technology Policy ISTECS Journal Vol. VIII/2006
43
ISTECS JOURNAL, Vol. VIII (2006) 44 – 53
Enhancement of Agricultural Product Quality Through Management on Controlled Environment M. Affan Fajar Falah Doctoral Student, United Graduate School Of Agricultural Sciences, Ehime University Department of Agroindustrial Technology, Fac. Agric.Technology UGM Jl.Sosio Yustisia Bulaksumur No: 1 Yogyakarta E-mail:
[email protected]
ABSTRACT Agricultural product quality can be improved through controlled environment during preharvest production and postharvest storage. This improvement is needed to reach costumer satisfaction with improve and keep the nutritional status of agricultural product. Enhanced the agricultural product quality must be support using some regulation by the government with several objectives, such as supplied good quality and safety of agricultural product that determined by research, securing food security, increase in welfare of the farmers and for strengthen the micro and the macroeconomic. Keywords : Agricultural product quality, Controlled environment, Costumer satisfaction, Government policy.
1. INTRODUCTION The definition of term Product Quality has been subject of large number of studies and some leader that contribute to the quality theory such as Edward Deming with the 14 points of quality management with his idea that known as `continual never ending improvement`, Joseph M Juran with Juran`s Trilogy (planning, quality and improvement) that defined the quality as fitness for use and also a Japanese leader for quality improvement in Japan, Kaoru Ishikawa, that also known as the inventor of the seven tools of quality control and his some philosophy in quality are `the quality begins with education and ends with education`, the first step of quality is to know the requirement of the consumer and the ideal state of quality control is when inspection is no longer necessary [9]. The other contributors for quality theory also important for comprehensive theory in quality based on his or her point of view that range from engineering, agriculture, statistics, management and also social science. On the other hand, quality of agricultural product can be defined with a specific term as a combination of attributes or properties that determined their value to the consumer with quality parameters that include intrinsic value such as appearance, texture, flavor, shelf life and nutritional value and depend on the whether the product is eaten fresh, processed or cooked and extrinsic factor that include production system for that product is producing, packaging, specific processing technology, fertilizer that used and also biotechnology to modifies properties of product [12]. Extrinsic factors Science and Technology Policy ISTECS Journal Vol. VIII/2006
44
can affect consumer acceptance of the product not internal or physical quality of the product. In developing countries, consumer demands now not only need agricultural product that look nice from appearance but also have an appropriate for texture, free of contaminants (for organic product), safety, contain sufficient of nutritional and health promoting substances [16]. Human needs for good quality of agricultural product for his or her health for everyday activities in their life. Agricultural product as a source of food that consume by the human is recommend the nutrition that must be eaten. The essential nutrients that used by human nutritionist include both of vitamins and essential nutrients, but plant nutritional scientist define only essential nutrients, because plants being autotrophs, in this paper we used term of nutritional elements based on human nutritionist. The human require at least 49 nutrients to meet their metabolic needs that shown in Table 1 [21]. Producing good quality of agricultural product can be enhanced through plant breeding with modification of genetic engineering. The classical approach for breeding cultivars is to select suitable phenotypes or mutants, which are then crossed, selfed, cloned or combined with populations depending on their reproductive biology. One of the known method is biofortification that define as increasing the bioavailable micronutrient contents of agricultural product through genetic variations through plant breeding. These tools is believed by the several scientists should be exploited by the nutritionist and public health communities to reduce malnutrition in the developing world that caused nutrition deficiency [21].However, in this paper we want to overview the other way to enhance agricultural product quality through management on controlled environment for preharvest production and postharvest storage. Table 1. The 49 known essential nutrients for sustaining human life No
Nutrients
Name of nutrients
1
Water and energy
water and carbohydrates
2
Protein
histidine, isoleucine, leucine, lysine, methionine, phenylalanine, threonine, tryptophan and valine
3
Lipids-fat
linoleic acid and linolenic acid
4
Macro-elements
Na, K, Ca, Mg, S, P, Cl
5
Micro-elements
Fe, Zn, Cu, Mn, I, F, B, Se, Mo, Ni, Cr, V, Si, As, Sn and Co
6
Vitamins
A, D, E, K, C (ascorbic acid), B1 (thiamin), B2 (ribovlavin), B3 (pantothenic acid), niacin, B6 (pyridoxal), folate, biotin and B12.
Source: [21]
Management on controlled environment can affect and improve quality of agricultural products. During pre-harvest production, control on environment condition such as temperature and humidity that influence the nutritional value of product, physical, mechanical and chemical composition are important. Location, season and also cultural practices affect the product quality of agriculture. For example, temperature can influence uptake and metabolism of mineral nutrients metabolism [1][12], daily low Science and Technology Policy ISTECS Journal Vol. VIII/2006
45
temperature can increase content of beta carotene in broccoli [19], low root temperature also increase in ascorbic acid content on welsh onion and spinach [10]. On the other hand, during post-harvest handling and storage, quality of agricultural product also strongly influenced by various environmental conditions, such as temperature, humidity and atmosphere conditions of the product. For example, low temperature can enhance storability and shelf life of most vegetable product; moisture loss from harvested agricultural product causes changes in structure, texture and appearance of the product; reduction of water loss from the product during postharvest period is critical control point for freshness and quality; and modified and controlled atmospheres are important for delay quality deterioration, delay vernalization and sprouting [6] [13] [16]. This overview will explained the processes and practices from preharvest production and postharvest storage for enhance and maintain quality of agricultural product. 2. QUALITY ENHANCEMENT DURING PREHARVEST PRODUCTION Controlled Environment Agriculture (CEA) is a combination of horticultural and engineering techniques that optimize crop production, crop quality, and production efficiency. CEA crops will probably be grown using some form of hydroponics, whether this be directly in liquid media, having the roots in contact with a thin film of liquid, having the roots in contact with a porous pipe that carries media, or in some porous matrix like vermiculite. Plants are maintained in controlled environment where environmental conditions, like lighting, nutrient supply, temperature and humidity, can be strictly controlled by computer. Well-managed, CEA operations can provide fresh produce (as well as flowers or pharmaceutical plants) of high quality and free of agriculture chemicals. CEA facilities can also be located in urbanized areas, thus not requiring the conversion of open or agricultural land. The two most important environmental variables for growing plants are temperature and light. Both parameters must be controlled to be uniform from one location to another in a greenhouse, and consistent from day to day. Many types of CEA facilities are common in use for agricultural scientist or for agricultural industry and farmers, such as greenhouse, phytotron, growth chamber and Controlled Ecological Life-Support System (CELSS). The most common and widely used of CEA is the greenhouse. The greenhouse is a defined as a transparent enclosure designed to grow or temporarily protect plants or building which should provide the conditions for satisfactory plant growth throughout the year. The growth factors such as light from solar radiation, temperature, humidity and air composition should be delivered and maintained at optimum levels for plant growth [7].The purpose of the good greenhouse, especially in construction, is to create the necessary climatic conditions for plant production throughout the year and to control the climatic factors as far as possible to the specified optimum. For these purposes, the greenhouse are required to have high light transmission, low heat consumption, sufficient ventilation efficiency, adequate structural strength, low cost construction and also low operating cost [22]. In plant production systems that grown in CEA , such as a greenhouse, , plants are frequently
exposed to the environmental stresses (high and low light stress, high and low temperature stress, water Science and Technology Policy ISTECS Journal Vol. VIII/2006
46
stress, low dissolved O2 stress, nutrient deficiency, etc). Environmental factors far from their optimum conditions bring the environmental stresses to the plants. That is, for the optimum cultivation in plant production systems, it is necessary to understand the physiological mechanisms of induction of environmental stresses and tolerance to the stresses in plants, and to establish the methods for removing and applying the environmental stresses in plant production systems [10]. Management of environmental stresses for plant food agriculture based on physiological basis can be shown in the figure 1 (Laboratory Bioproduction System, Kochi University). Environmental Stresses
Plant Production Systems
Effects to plants
Methods and Equipment
• High Light Stress Removing
• Low Light Stress
• Removing of High and Low Light Stress − High Yield, Labor Saving
• Double Seesaw Mechanics Production
• Strengthening of Root Physiological Function
• O2 Application to Solution • Storing Heat to Under Ground
• High Temperature Stress
• Treatment of Short Term Stress • Production, Storing, Acclimation of Seedling
• Low Temperature Stress • Water Stress
Physiological responses to stress Application
• Induction of Osmoregulation
• Salt Stress
• Induction of Antioxidation
• Low DO2 Stress
• Induction of Stress Tolerance
Evaluation of Physiological Effect
• Removing of High and Low Temperature stress with Energy Saving • Accumulation of Sugar and Amino Acid • Accumulation of Antioxidants • Hardening of Seedling
Figure 1. Schematic diagram of the removing and applying the environmental stress in plant production systems.
We applied the high temperature stress on root by using a Nutrient Film Technique in the greenhouse with tomato plants. We found that water and nutrients absorption firstly increased for short-term effect of high temperature on root, just like hysteresis effect, and the long term effect of high temperature on root of water and nutrients absorption by root and translocation to the part of plant (roots, leaves, stem and fruits) during harvested for tomato production were depressed and lower content of nutrient compare than those in the plants that grown in the control temperature [1] [2] [3]. However, when we applied the high temperature on root during salt stress condition, we did not found the differences of nutrients content in the part of organ compare than that grown in control temperature with salt stress [3]. On the other hand, we also applied salt stress using a deep seawater for improve production of tomato plants. Tomato plants (Lycopersicon esculentum Mill.) were grown in the NFT system, where the concentrated deep seawater was applied for the short-term salt stress treatment for only two weeks at the stage of rapid fruit growth. From these physiological analyses, it verified that the short-term application of the concentrated deep seawater at the stage of rapid fruit growth can induce the osmoregulation in the phloem transport to fruits and can produce high quality tomatoes enriched in sugar, acid, minerals and good flavor without occurrence of extremely small-sized fruits and blossom-end rot [20]. For low temperature on roots, we also applied to the welsh onion and spinach during summer season Science and Technology Policy ISTECS Journal Vol. VIII/2006
47
in the greenhouse with Deep Flow Technique. In the welsh onion and spinach, low temperature were inhibited water and nutrient transport from the root, however dry matter content of nutrients were higher compare than control temperature on root [10]. On the other hand, when we applied the low light stress of solar radiation in the strawberry cultivation through the double seesaw mechanic system, we found that productivity and appearance of the strawberry fruit were increase and look more fresh [11]. Thus, using the many kinds of plant instrumentations enabled the assessment of environmental stress in plant production systems. In future, it is expected that the establishment of the other methods for removing and applying the environmental stress on the basis of results of the present study, for example removing the high root temperature stress by using underground system with energy saving from solar cell in the greenhouse. Furthermore, the establishment of the methods for removing and applying the environmental stress will contribute to the achievement of the optimum cultivation in plant production systems for high quality product of agriculture. 3. QUALITY ENHANCEMENT DURING POSTHARVEST STORAGE The six primary environmental variables usually controlled in storage of postharvest product to maintain quality are duration (shelf life) , temperature, relative humidity, and the concentrations of O2, CO2, and ethylene. The ‘optimum’ storage environment for each commodity is designed to maintain these variables within a set of limits that produces the maximum storage life for most of the individual members of the commodity. These storage variables could be thought of as a volume in six-dimensions. Agricultural product of the major cultivars and the major growing locations have specific recommendations that often differ significantly from one another. Yet even within these refined storage conditions, agricultural product from different growing seasons and parts of the tree, and those experiencing different cultural practices and harvesting techniques can respond differently to the ‘optimum’ storage environment. Other environmental variability that influences the response of the commodity during storage include microbial load, light, orientation of the product in the gravitational field, and the concentration of other gases; the most important being ethylene [18]. The history of storage using modified atmospheres is long storied and ultimately, extremely successful. These technologies can be segregated into two classes based on the manner in which the atmospheres ore obtained and maintained. One is referred to as modified atmosphere packaging (MAP), in which a barrier passively limits gas exchange and the living produce enclosed in the package containing the barrier impacts the headspace atmosphere by virtue of its respiratory activity. In MAP, for instance, the atmosphere can be initially modified, the barrier may be gas permeable films or perforated films, or both, gases may be actively released or scavenged in the package, a partial vacuum can be imposed, the package may be used at the distributor level only, the consumer level only, or both, there may be sensors incorporated in the package, and so on [6] [13]. The second is termed controlled atmosphere (CA) storage, in which the atmosphere is either manually or mechanically controlled to achieve target headspace gas concentrations. In CA storage, there are also numerous variations in design and application, ranging from atmosphere establishment by the respiration of the product to active modification by flushing, scrubbing, or vacuum (in low pressure Science and Technology Policy ISTECS Journal Vol. VIII/2006
48
storage). Sense-and-respond systems have been designed and installed that use the responses of the stored product to assist in decision-making regarding atmosphere application. Nevertheless, the aim of MAP and CA technologies is to take advantage of physiological responses of the enclosed plant material and/or accompanying pathogens to modification of the atmospheric content of respiratory gases O2 and CO2. The physiological processes typically targeted for control via atmosphere modification and the impact of low O2 on ethylene biology is probably the most economically important application of atmosphere modification and is commonly applied using CA facilities [6] [13]. The air composition can be controlled using CA and MAP facilities, because of the availability of lower O2 and high CO2 are additive in the storage atmosphere and vary between the agricultural commodity. Low concentration of O2 reduces respiration, ethylene synthesis and action, but it also stimulates anaerobic respiration, the production of off flavors, and possible microbe growth. Furthermore, increased CO2 concentration are difficult to explain, but normally stimulates respiration in stored roots and bulb, inhibit enzymes in Kreb`s cycle and decreasing pH and the cell sap [16] [17]. A number of volatile inhibitors of ethylene action (e.g. 1-MCP, NBD, PPOH) have been discovered in the past decade. Pre-storage treatment with 1-MCP can render the commodity insensitive to ethylene. Inhibiting ethylene action can have tremendous commercial benefits for the storage of sensitive commodities. However, not all ripening changes are controlled by ethylene and some ethylene responses may retard disease development. Additionally, ethylene production by some non-climacteric commodities (e.g. citrus, vegetative tissue) is under negative feed-back control, so that reducing perception of ethylene will actually stimulate its production. Also, the storage life of many ethylene sensitive commodities is not terminated by their response to ethylene, but by water loss, mechanical injury, or disease development. Use of compounds like 1-MCP will not compensate for poor temperature control, handling practices, or sanitation [18]. Temperature is the all-important consideration for postharvest food quality, and the principle of overriding importance is to ensure a proper environment with the correct temperature range in food storage. Stores with non-uniform air and temperature distribution often suffer from condensation, consequent rotting and poor quality of fresh foods . A uniform temperature distribution, on the other hand, can lead to minimum deviations. Lower temperatures than required can cause chilling injuries and higher temperatures accumulated in certain points of the storage area will result in greater postharvest quality loss and increased respiratory and microbial activity and hence shortened shelf life [16]. Low temperature can enhance the storability and shelf life of most fresh agriculture product as long as the temperature is above the threshold level of chilling injury or freezing damages of the tissue. With the decreasing temperature, the kinetic energy of water molecules diminishes, resulting in a reduced vaporization rate of liquid water with a lower transpiration and water loss [18]. On the other hand, prestorage heat treatments to control decay are often applied for a relatively short time (minutes), because the target pathogens are found on the surface or in the first few cell layers under the skin of the fruit or vegetable. Heat may be applied to fruit and vegetables in several ways: by hot water dips, vapor heat, or hot dry air or by hot water rinsing and brushing. Vapor heat treatment was developed mainly for insect control, while hot dry air has been used for both fungal and insect control. Water is the preferred Science and Technology Policy ISTECS Journal Vol. VIII/2006
49
medium for most applications since it is a more efficient heat transfer medium than air. Heat treatments can also be used to inhibit ripening processes or to induce resistance to chilling injury and external skin damage during storage, thus extending storability and marketing, and treating with hot water has become increasingly accepted commercially, and a significant improvement has been made with the addition of brushing [8]. 4. RECENT POLICY FOR DEVELOP PRODUCT QUALITY IN AGRICULTURE From Indonesian Agriculture Development Plan 2005-2009 [4] explained that in general, position of Indonesian technology status for several agriculture commodities are relatively left behind as compared to other ASEAN countries, for estate crops commodities Indonesian technology is left behind as compared to Malaysia, while horticulture is left behind as compared to Thailand. For food processed product, Indonesian product is relatively left behind as compared to Thailand and Vietnam. Those are influenced by more consistent of government attention in developing agribusiness channels from up-stream, down-stream up to marketing system for both processed and fresh products. The vision of agricultural development for the period of 2005-2009 is to realize strong agriculture for strengthening food security, improvement of value added and competitiveness of agricultural products and the improvement of farmer welfare. Strong agriculture is agriculture sector which is able to grow sustainability and having characteristics as follows: (1) knowledge is a main foundation in decision making, strengthening intuition, custom and tradition; (2) advance technology is a main instrument in resource utilization; (3) market mechanism is a main media in transaction of goods and services; (4) efficiency and productivity are basic of resource allocation; (5) quality and competitiveness is an orientation, expression and objective; (6) professionalism is a main character; (7) engineering is center of added value, so that every produced product always meet determined requirements. There are three main targets of agricultural development which must be reached in the next five years namely: (1) the improvement of national food security covering improvement of production capacity of agricultural commodities and decreasing the dependency to food import around 5-10 percent of domestic demand; (2) the improvement of value added and competitiveness advantage of agricultural commodities covering the improvement of the qualities of agricultural products, the improvement agricultural product processing diversification and the increase export and export surplus of agricultural product; and (3) the improvement of farmer welfares. Policy on improving promotion and protection on agricultural commodities is directed to: (a) Formulating policy on subsidy, input production, output prices, interest rates and credit for agriculture activities; (b) Increasing export and restricting import; (c) Policy on import tariff and import regulation; (d) Increasing productivity and efficiency an agricultural activities; (e) Qualities improvement and product standardization through implementing production technology, postharvest management; and (f) Empowerment of marketing system and protecting agricultural activities. We can find that the real condition of the Indonesian agricultural product, left behind the ASEAN countries. The government through the ministry of agriculture have some policy to increase the quality of agricultural product for strengthen the micro and macroeconomic and to improve welfare of Science and Technology Policy ISTECS Journal Vol. VIII/2006
50
the farmers. We can learn from any countries to build the strong agriculture product, for example we can learn from Japan for keep the quality of agricultural product and improve the welfare of the farmers. Although, Japan food self-sufficiency on calori basis is only 40%, but for production basis is 70% and the quality of their product is one of the good quality product that can be produced through good quality production in the field and also greenhouse production. The assuring of their quality product was support by the policy from the governmental office (i.e ministry of agriculture, fisheries and forestry) through some policy in Japan Agricultural Standard that was defined, firstly from customer voice and strengthen by the policy on the research and technology that support this program. Furthermore, their policy to improve the product quality, although rise in price, have the main purpose to increase the welfare of the farmers and stabilize the micro and macroeconomic. One of the famous program launched to promote Japan high quality of agricultural product is the Oishii (Delicious) campaign that start from March 2005 [5]. The government said that Japan brand of agricultural product may command high prices, but they guarantee and assure the finest quality and the very best flavors, it indeed the real taste of Japan. This declaration is very important for promoting and marketing of their product. 5. CONCLUSION Human needs for good quality of agricultural product for his or her health for everyday activities in their life. Agricultural product as a source of food that consume by the human is recommend the nutrition that must be eaten. This nutrition can be improved through some ways or methods, one of the method is controlled environment during preharvest production and postharvest storage. By research, this method can made clearly explained the mechanism and the regulation that need to enhance and maintain the nutritional status of plant product that human needs. Government regulation and policy is needed to support this program, not only for improve and maintain the agriculture product quality but also for the other reason such as securing food security, increase in welfare of the farmers and for strengthen the micro and the macroeconomic. To assure the quality of agricultural product, the simultaneous ways from research and government policy to implement the quality program for farmers are important to do. We hope the delicious agricultural product from Indonesia can produce with the high quality and safety. AKNOWLEDGEMENT Authors depply grateful thanks to Prof Masaharu Kitano for his guidance, and for members of Laboratory Bioproduction System in Kochi University, especially for Wajima Takahiro and Hidaka Kota that gave their data. REFERENCES [1] Affan, M, F, F., High Temperature Effect on Root Absorption in Hydroponic System, Master Thesis, Kochi University (2004), pp 90. [2] Affan, M, F, F., Kitano, M., Yasutake, D., Wajima, T. Analyses of High Temperature and Salt Stress Science and Technology Policy ISTECS Journal Vol. VIII/2006
51
on Roots in Tomato Hydroponic, In Proceeding of International Conference on Research Highlight and Vanguard Technology on Environmental Engineering in Agriculture System. Kanazawa (2005a), 215-220. [3] Affan, M, F, F., Kitano, M., Yasutake, D., Wajima, T., Yasunaga, T. Study on Root Absorption responding to Environmental Stress by Using Hydroponic Systems. Phyton Annales Rei Botanica (2005b), 223-228. [4] Anonim, Indonesian Agriculture Development Plan 2005-2025. Ministry of Agriculture Republic of Indonesia (2005). [5] Anonim, Basic Plan and Overview. Ministry of Agriculture , Fisheries and Forestry Japan. (2006). [6] Beaudry, R., Luckanatinvong, V., Solomos, T. Maintaining Quality with CA and MAP. Acta Horticulturae 712 ( 2006), 245-252. [7] Enoch, H, Z., Enoch,Y. The History and Geography of The Greenhouse. In : Greenhouse Ecosystem, eds : Stanhill, G and Enoch, H, Z., Amsterdam : Elsevier (1998), 1-16. [8] Fallik, E. Prestorage Hot Water Treatments (Immersion, Rinsing and Brushing). Postharvest Biology and Technology 32 (2004) 125–134. [9] Foster, S,T., Managing Quality : An Integrated Approach, New Jersey: Pearson Prentice Hall, (2004). [10] Hidaka, K., Assesment of Environmental Stress in Plant Production System. Master Thesis of Kochi University (2006). Pp 109. [11] Hidaka, K., Ito, E., Kitano, M., Imai, S., Innovative Cultivation of Strawberry on Vertically Moving Beds Controlled by Double Seesaw System. In Proceeding of International Conference on Research Highlight and Vanguard Technology on Environmental Engineering in Agriculture System. Kanazawa (2005), 257-262. [12] Kader, A, A., Quality Parameters Of Fresh-cut And Vegetable Products, In Fresh-cut Fruits and Vegetables: Science, Technology and Market Ed: Lamikanra, O., Florida : CRC Press (2002), pp. 11-20. [13] Kader, A, A., Rolle, R, S. The Role of Postharvest Management in Assuring The Quality and Safety of Horticultural Produce, Rome : FAO (2004), pp 52. [14] Jongen, W, M, F., Food Supply Chains : From Productivity Toward Quality, In Fruit and Vegetable Quality : An Integrated View, Eds : Shewfelt, R, L. and Bruckner, B., Pensylvania : Technomic Publishing Company Inc (2000), pp 3-18. [15] Nilsson, T. Postharvest Handling and Storage of Vegetables.In Fruit and Vegetable Quality : An Integrated View, Eds : Shewfelt, R, L. and Bruckner, B., Pensylvania : Technomic Publishing Company Inc (2000), pp 96-116. [16] Nicolai, B, M., Beullens, K., Bobelyn, E., Hertog, M, L, Schenk, A., Vermeir, S, Lammerty, J. Systems to Characterise Internal Quality of Fruit and Vegetables. Acta Horticulturae 712 ( 2006), 59-65. [17] Ozcan, S, E., Cangar, O., Vranken, E., Berckmans, D. Predicting 3D Spatial Temperature Uniformity in Food Storage Systems from Inlet Temperature Distribution. Postharvest Biology and Science and Technology Policy ISTECS Journal Vol. VIII/2006
52
Technology 37 (2005) 186–194 . [18] Saltveit, M. Is it Possible to find an optimal controlled atmosphere?, Postharvest Biology and Technology 27 (2003) 3 -/13. [19] Schreiner, M., Huskens-Keil, S., Krumben, A., Schonhof, I., Linke, M. Environmental effects On Product Quality.In Fruit and Vegetable Quality : An Integrated View, Eds : Shewfelt, R, L. and Bruckner, B., Pensylvania : Technomic Publishing Company Inc (2000), pp 85-94. [20] Wajima, T., Kitano, M., Araki, T., Jun, Y., Ishikawa, K., Matsuoka, T. High Quality Tomato Production by Suitable Application of Concentrated Deep Seawater. In Proceeding of International Conference on Research Highlight and Vanguard Technology on Environmental Engineering in Agriculture System. Kanazawa (2005), 83-88. [21] Welch, R, M., Graham, R, D., Breeding for Micronutrients in Staple Food Crops from A Human Nutrition Perspective, Journal Of Experimental Botany 55 (2004), pp 353-364. [22] Zabeltitz, C, V. Greenhouse Structure, In : Greenhouse Ecosystem, eds : Stanhill, G and Enoch, H, Z., Amsterdam : Elsevier (1998), 17-70.
Profile M.Affan Fajar Falah was born in Kendal, Central Java Indonesia on April 10th, 1975. He is a staff of Gadjah Mada University, Faculty of Agricultural Technology, Department of Agro-industrial Technology, Jogjakarta from 1999. Recently, he has been pursuing a PhD program at Ehime University Faculty of Agricultural Sciences, Japan. He completed a B.Sc degree in 1998 from same Department at Gadjah Mada University, Jogjakarta Indonesia and Master`s Degree on Agriculture in 2004 from Kochi University, Japan. His current interest is in Controlled environment agriculture during production and storage and quality improvement.
Science and Technology Policy ISTECS Journal Vol. VIII/2006
53
