Advanced
Design of Nonlinear Observer for Flexible Joint Manipulator with Nonlinear Stiffness Based on Acceleration
Design of Nonlinear Observer for Flexible Joint Manipulator with Nonlinear Stiffness Based on Acceleration
The Transactions of The Korean Institute of Electrical Engineers. 2015. Mar, 64(3): 451-457
Copyright © 2015, The Korean Institute of Electrical Engineers
This is an Open-Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/3.0/)which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.
  • Received : November 28, 2014
  • Accepted : February 27, 2015
  • Published : March 01, 2015
Download
PDF
e-PUB
PubReader
PPT
Export by style
Article
Author
Metrics
Cited by
TagCloud
About the Authors
승 준 이
Korean Institute of Industrial Technology, Korea.
형 종 김
ASRI, Dept. of Electrical and Computer Engineering, Seoul National University, Korea, and also Korean Institute of Industrial Technology, Korea.
경 태 남
Corresponding Author : Korean Institute of Industrial Technology, Korea. E-mail:robotnam@kitech.re.kr
태 용 국
Dept. of Electronic Eng., Sungkyunkwan Univerity, Korea.

Abstract
In this paper, we consider the observer design problem that truly reflects the nonlinear stiffness of the manipulators. The two key ideas of our design are that (a) estimation error dynamics of the manipulator equipped with accelerometer dose not dependent on nonlinearities at the link part, when the measured signals are the motor position and the output of the accelerometer and (b) the nonlinear stiffness is indeed a Lipschitz function. In order to effectively compensate the nonlinear stiffness, the gain of the proposed observer is carefully chosen from the ARE(algebraic Riccati equations) which depend on Lipschitz constant. Comparative simulation result verifies the effectiveness of the proposed solution.
Keywords
1. 서 론
유연관절로봇(flexible joint manipulator)은 자동화 공정을 통하여 효과적인 생산성을 필요로 하는 산업현장에서 널리 쓰이고 있다 [1] [4] . 따라서 정밀하고 신속한 움직임을 얻기 위해 시스템에 적용되는 제어기의 중요성이 증대되고 있다. 이런 유연관절로봇의 정밀한 위치제어 문제는 오랜 기간 연구되어왔으며, 많은 효과적인 제어방법들이 제안되었다 [1] - [6] . 하지만, 대다수의 제어기들이 로봇의 정확한 상태(위치, 속도 등)정보를 필요로 하는 것에 반해, 로봇 상태 추정에 대한 연구는 많이 부족한 실정이다. 또한 몇몇 제안된 관측기(observer) 역시 유연관절로봇의 비선형 특성이나 모델 불확실성으로 인하여 정확한 로봇상태를 알아내기가 쉽지 않다. 특히, 유연관절의 특성상 모터의 상태만 관측 가능한 제약 속에서 링크단의 상태를 알아내는 일은 쉬운 일이 아니다 [3] [10] .
최근에 보다 정확한 로봇 상태추정을 위하여 가속도(acceleration) 정보를 기반으로 한 관측기가 제안되었으며, 여러 연구가 진행되고 있다 [11] - [15] . 이는 사용 및 운용이 쉽지 않은 위치 센서 대신 로봇 링크(link) 쪽에 가속도 센서를 장착하여 링크 쪽 비선형특성의 정보를 얻어 관측기를 설계함으로써, 효과적으로 로봇의 상태를 추정한다. 일반적으로 이런 연구에서는 로봇의 선형강성(linear stiffness)만을 고려하여 설계되었다. 하지만, 실제 로봇의 강성은 비선형특성을 가지고 있어 선형강성만을 고려한 기존 관측기만으로는 정밀한 상태추정이 용이하지 않다.
본 논문에서는 로봇의 비선형강성(nonlinear stiffness)을 고려한 가속도기반의 비선형 관측기를 제안한다. 제한한 관측기는 비선형강성을 가진 로봇 모델을 Lipschitz 비선형시스템으로 변형하여 설계되었으며, 그 관측기의 이득(gain)은 관측오차를 점진적으로 0으로 수렴할 수 있도록 설계된다. 끝으로 실제 로봇모델을 이용한 모의실험을 통하여 기존의 선형강성만을 고려한 관측기보다 향상된 상태추종 성능을 확인한다.
2. 본 론
먼저 유연관절로봇의 비선형강성 특성에 대해서 소개를 하고, 2.2절에서 가속도정보를 사용한 비선형 관측기의 제시와 안정도 분석을 한다. 그리고 2.3절에서 관측기 이득 설계 시 도움이 되는 Lipschitz 상수를 줄이는 좌표변환(coordinate transformation)을 소개한다.
- 2.1 비선형강성
유연관절로봇의 특성 상 각축간의 기어(gear)나 기구적인 특성으로 인한 강성(stiffness)이 존재하며, 이런 강성은 실제로 아래 그림 1 과 같은 토크(torque)와 링크단과 모터단의 비틀림 각도 ( θ )의 비선형함수의 특성을 가지고 있다 [7] - [9] .
PPT Slide
Lager Image
비선형강성 Fig. 1 Nonlinear stiffness
이런 비선형함수를 포함하는 비선형강성은 다음과 같은 식으로 표현될 수 있다.
PPT Slide
Lager Image
여기서 θ 는 비틀림 각도며, k 1 k 2 는 선형과 비선형 상수이다. 또한 θB 는 굴절중단점(breakpoint deflection)이다. 우리는 이런 비선형강성을 갖는 유연관절로봇을 위한 상태관측기를 설계한다.
- 2.2 관측기 설계
본 논문에서는 이관성(two-inertial) 시스템의 유연관절로봇 모델을 기반으로 그 시스템의 상태를 관측하기 위한 관측기를 설계한다. 즉, 다음과 같은 미분방정식으로 이관성시스템을 표현할 수 있다 [7] .
PPT Slide
Lager Image
여기서 θl ωl 은 링크의 위치와 속도를 나타내며, θm ωm 모터의 위치와 속도를 나타낸다. Jl ( θl )은 링크관성, D 는 댐핑, C ( θl , ωl )는 코리올리와 원심력, G ( θl )는 중력, Jm 은 모터관성, 그리고 τ 는 입력토크를 나타낸다. K ( θl , θm )는 비선형강성을 나타내는 것으로 다음과 같이 식(1)의 비틀림 각도( θ )를 θm θl 로 바꾸어 표현될 수 있다.
PPT Slide
Lager Image
여기서 k 1 k 2 는 선형강성과 비선형강성을 각각 나타내며, θB 는 굴절중단점이다. θd : = θm θl 을 정의하면, 위 식은 다음과 같이 적을 수 있다.
PPT Slide
Lager Image
여기서
PPT Slide
Lager Image
PPT Slide
Lager Image
이다. 또한 y 2 는 링크의 가속도를 나타내는 것으로 우리는 다음과 같은 가정을 한다.
가정1: 링크 각가속도를 나타내는 식(4)의 y 2 는 측정가능(measurable)하다. ◇
위와 같은 시스템으로부터 우리는 아래와 같은 사실을 알 수 있다.
보조정리1: 식(4)의 함수 𝜙 3 ( x )는 globally Lipschitz이다. 즉, 다음 식을 만족하는 상수 γx > 0이 존재한다.
PPT Slide
Lager Image
증명: 𝜙(•)는 연속이고 미분 가능한 함수이다. 또한, φ 1 ( x 3 ), φ 2 ( x 3 ), 그리고 φ 3 ( x 3 ) 는 각각 연속인(continuous) 함수이며, φ 1 (− θB ) = φ 2 (− θB )와 φ 2 (− θB ) = φ 3 (− θB )이다. 그리고
PPT Slide
Lager Image
는 다음과 같다.
PPT Slide
Lager Image
여기서
PPT Slide
Lager Image
,
PPT Slide
Lager Image
, 그리고
PPT Slide
Lager Image
역시 각각 연속인 함수이며,
PPT Slide
Lager Image
PPT Slide
Lager Image
이다. 이로부터 아래의 식과 같이
PPT Slide
Lager Image
는 고르게 유계(uniformly bounded)이다.
PPT Slide
Lager Image
위와 같이 함수 φ 3 ( x )는 연속이며 고르게 유계이기 때문에
PPT Slide
Lager Image
을 만족하는 어떤 γx > 0가 존재한다. 그러면,
PPT Slide
Lager Image
이며 이로부터
PPT Slide
Lager Image
을 알 수 있다. 따라서 함수 𝜙(•)는 globally Lipschitz이다. (좀 더 자세한 내용은 [18] 의 보조정리 3.3에서 확인 할 수 있다.) ■
이제,
PPT Slide
Lager Image
의 Lipschitz 상수 γx 를 다음과 같이 3 k 2 θ 2 B 로 얻을 수 있다. 첫째, x 3 ≺ − θB 구간에서 다음과 같은 부등식을 얻을 수 있다.
PPT Slide
Lager Image
다음으로 | x 3 |≤ − θB 구간에서는 다음과 같으며,
PPT Slide
Lager Image
마지막으로 x 3 ≻ − θB 구간에서는 다음을 얻을 수 있다.
PPT Slide
Lager Image
주목1: 식(4)와 보조정리1로부터 우리는 다음과 같은 사실을 알 수 있다. 만약 행렬 A 의 인자로 비선형 강성 K ( θl , θm )의 k 1 을 사용하지 않으면 행렬 짝 ( A , C )는 관측가능하지 않으며(unobservable), 그 경우 일반적으로 관측기 설계가 용이하지 못하다 [19] . 따라서 본 논문에서는 시스템 (2)을 (4)로의 변환을 통하여 강성 K ( θl , θm )의 k 1 부분을 이용하여 시스템을 관측가능(observable)하게 만들었으며, 또한 나머지 강성부분
PPT Slide
Lager Image
를 Lipschitz 조건을 만족하도록 만들었다. ◇
지금까지의 결과로 시스템(4)는 다음과 같은 Lipschitz조건을 만족하는 비선형시스템으로 생각할 수 있다.
PPT Slide
Lager Image
여기서
PPT Slide
Lager Image
이다.
이제, 위의 비선형강성을 갖는 유연관절로봇에 대한 상태추정을 위하여 다음과 같은 관측기를 제안한다.
PPT Slide
Lager Image
여기서 L 은 관측기 이득(gain)을 나타내는 것으로 비선형함수
PPT Slide
Lager Image
을 위하여 그 설계방법을 아래에서 제시 할 것이다. 우선, 관측성능을 살펴보기 위하여 관측오차
PPT Slide
Lager Image
을 정의하면, 그것의 오차(error)시스템을 다음과 같이 구할 수 있다.
PPT Slide
Lager Image
여기서 Aob = A LC 이다.
식(13)에서 비선형부분
PPT Slide
Lager Image
이 없는 일반적인 선형관측기 ( ė = ( A LC ) e )에서는 그 수렴성을 보이기 위하여 관측기 이득 L 을 행렬 A LC 이 Hurwitz하게만 설계하면 된다 [19] . 하지만 우리는 식(13)와같이 비선형부분이 존재하기 때문에 관측기 이득 설계에 좀 더 다른 방식을 필요로 하게 된다. 우선, 보조정리1로부터 위 식(13)의 비선형부분이 Lipschitz 조건을 만족하는 것을 알 수 있으며, 이를 위한 관측기 이득 설계 방법에 대한 연구는 기존에 다양하게 진행되었다 [16] [17] . 본 논문은 ARE(Algebraic Riccati equations)방식 [16] 을 이용하여 관측기 이득을 설계한 방법을 소개한다.
우선, 어느 ε ≻ 0에 대하여 양의 정부호행렬(positive definite matrix)인 해(solution) P 를 갖는 다음과 같은 ARE를 생각하자.
PPT Slide
Lager Image
위 식으로부터 관측기 이득은 다음과 같이 설계 가능하다.
PPT Slide
Lager Image
위와 같이 설계된 관측기와 관측기 이득에 대하여 다음과 같은 주요 정리를 제시한다.
정리1: 주어진 유연관절로봇 (2)과 설계된 관측기 (12)에 대하여, 만약 ARE 식(14)에 양의 정부호행렬 해 P 을 가지고 관측기 이득을 식(15)과 같이 설계한다면 관측기 오차 (13)는 점진적으로(asymptotically) 0으로 수렴한다. ◇
증명: 우선, 식(14)와 (15)로부터 다음을 알 수 있다.
PPT Slide
Lager Image
따라서 다음과 같은 부등식이 만족한다.
PPT Slide
Lager Image
그러면 [16] 의 보조정리1에 의해서 다음과 같은 부등식을 만족하는 어떤 양의 정부호행렬 P 1 이 있다는 것을 알 수 있다.
PPT Slide
Lager Image
이제, 시스템 (13)에 대한 안정도를 보이기 위해서 다음과 같은 Lyapunov함수후보를 생각한다.
PPT Slide
Lager Image
여기서 행렬 P 1 는 양의 정부호이다. 그러면 그 함수의 미분은 다음과 같이 구할 수 있다.
PPT Slide
Lager Image
따라서 식 (18)로부터
PPT Slide
Lager Image
을 알 수 있으며, e = 0은 시스템(13)의 평형점(equilibrium point)이다. 이제, 임의의 a > 0에 대한 집합 Ba = { e R | ∥ e ∥≤ a }를 생각하면, b = { e Ba V ( e ) ≤ b }⊂ Ba 를 만족하는 b > 0 를 선택할 수 있다. 식(20)로부터 집합 b 안에서 시작한 값은 시간이 지나도 집합 b 에 머물러 있는 것을 알 수 있다. 여기서 a 는 임의의 값으로 설정했기 때문에 모든 시간 t > T 에 대하여 ∥ e ( t )∥< a 인 어떤 T > 0 를 생각할 수 있으며, 이로 인해
PPT Slide
Lager Image
을 알 수 있다. (좀 더 자세한 내용은 [18] 의 정리 4.1을 참고 할 수 있다.) ■
- 2.3 Lipschitz 상수를 위한 좌표변환
식(14)의 해 P 는 Lipschitz 상수 γ 에 의존하는 것으로 만약 γ 가 아주 큰 값을 가질 경우 식(14)을 만족하는 해를 찾기가 쉽지 않다. 이는 ARE인 식(14)의 Hamiltonian 행렬의 고유치(eigenvalue)의 실수 값이 0에 가깝기 때문에 발생한다 [20] . 따라서 본 장에서는 [17] 에서 사용된 특정 좌표변환을 통하여 그 Lipschitz 상수 γ 값을 줄일 수 있는 방법을 소개한다.
우선, 임의의 β 값과 함께한 다음과 같은 변환 행렬을 생각하자.
PPT Slide
Lager Image
그러면, 변환 z := Tx 를 통해서 다음과 같은 시스템을 얻을 수 있다.
PPT Slide
Lager Image
여기서 각각의 기호는 다음과 같다.
PPT Slide
Lager Image
그리고 설계된 관측기(12)에 대해서도 변환
PPT Slide
Lager Image
을 적용하면 다음과 같은 시스템을 얻을 수 있다.
PPT Slide
Lager Image
그러면 새로운 오차
PPT Slide
Lager Image
에 대한 미분방정식을 다음과 같이 구할 수 있다.
PPT Slide
Lager Image
여기서
PPT Slide
Lager Image
이다. 그러면, 변환된 좌표에서 Lipschitz 상수 역시 다음과 같이 구할 수 있다.
PPT Slide
Lager Image
여기서 변환된 좌표에서 Lipschitz 상수는
PPT Slide
Lager Image
이다. 따라서 기존의 식(14)대신 다음과 같은 ARE를 만족하는 해
PPT Slide
Lager Image
를 구한다.
PPT Slide
Lager Image
여기서
PPT Slide
Lager Image
이다. 그러면, 그 해
PPT Slide
Lager Image
를 바탕으로 다음과 같이 관측기 이득을 설계할 수 있다.
PPT Slide
Lager Image
결국, 0보다 크고 1보다 작은 β 값을 이용하여 Lipschitz 상수값을 줄일 수 있다.
- 2.4 모의실험
본 논문에서 제안된 관측기의 효율성을 검증하기 위하여 Matlab / Simulink로 모의실험을 진행한다. 실제에 가까운 비교를 위하여 모의실험에 사용된 데이터와 파라미터는 [7] 에서 제시된 값을 사용하며, 대상 시스템의 주요 파라미터는 아래 표 1 과 같다.
대상 시스템 파라미터Table 1System parameters
PPT Slide
Lager Image
대상 시스템 파라미터 Table 1 System parameters
본 모의실험에서는 제시한 비선형관측기와 논문 [11] 의 선형강성만을 고려한 관측기와의 성능을 비교한다. 이를 위하여 표 1 과 같이 동일한 유연관절로봇의 파리미터를 사용했으며, 동등한 비교를 위하여 설계된 관측기 이득 역시 같은 값을 사용하였다. 아래 표 2 는 설계된 관측기의 설계 파라미터 값을 나타낸다. 여기서 Lipschitz 상수를 줄이기 위한 좌표변환의 β 값은 식(26)에 대응하는 Hamiltonian 행렬의 고유치의 실수 값을 0에 가깝지 않게 설계했으며, 표3 과 같이 β 값에 따라 그 고유치의 변화를 확인할 수 있다.
관측기 설계 파라미터Table 2Design parameters
PPT Slide
Lager Image
관측기 설계 파라미터 Table 2 Design parameters
Hamiltonian 행렬의 고유치Table 3Eigenvalue of Hamiltonian matrix
PPT Slide
Lager Image
Hamiltonian 행렬의 고유치 Table 3 Eigenvalue of Hamiltonian matrix
그림 2 는 실제 링크의 위치를 추종하는 결과를 보여준다. 확대된 오른쪽 그림에서 본 논문에서 제시한 관측기가 기존의 관측기보다 실제 링크 위치를 좀 더 정밀하게 추정하는 것을 볼 수 있다.
PPT Slide
Lager Image
링크 위치 추정 Fig. 2 Estimation of link angle
그림 3 은 링크위치( θl )에 대한 실제 측정값과 추정값의 오차를 나타낸 것으로 비선형강성을 고려한 관측기가 기존의 선형강성의 관측기보다 향상된 상태추종 오차 성능을 보여준다. 그리고 그림 4 는 비틀림각( θd = θm θl )에 대한 성능비교를 나타낸다. 이 결과 역시 본 논문에서 제시한 관측기가 향상된 상태추종 오차를 보여준다. 이와 같은 결과는 실제 유연관절로봇이 비선형강성을 포함하고 있기 때문에 이를 고려한 본 논문의 관측기가 기존의 선형강성만을 고려한 관측기보다 더 정밀하게 상태를 관측한다고 볼 수 있다.
PPT Slide
Lager Image
링크위치(θl) 추종 오차 Fig. 3 Estimation error of link angle θl
PPT Slide
Lager Image
비틀림 각(θd = θmθl) 추종 오차 Fig. 4 Estimation error of difference between link and motor angle θd = θmθl
3. 결 론
본 논문에서는 가속도정보를 이용한 유연관절로봇을 위한 상태관측기를 제시하였다. 제시한 관측기는 Lipschitz 비선형시스템 기반으로 설계되었으며, 그 안정도와 성능을 분석하였다. 특히, 본 연구는 기존의 접근과 달리 강성의 비선형성을 고려하여 실제 시스템에 좀 더 가까운 상태를 추정할 수 있었다. 따라서 일반적인 선형강성만을 고려한 관측기에 비하여 관측성능이 향상되었고, 이를 모의실험을 통해 확인 할 수 있었다.
Acknowledgements
본 연구는 지식경제 기술혁신사업 산업융합원천기술개발사업 (과제번호: 10041827) 지원에 의하여 이루어진 연구로서, 관계부처에 감사 드립니다
BIO
이 승 준(Seung-Joon Lee)
2000년 영남대학교 전기전자공학부 졸업.
2003년 동 대학원 전기공학과 석사졸업.
2003년∼2010년 ㈜LS산전 선임연구원.
2010년∼현재 한국생산기술연구원 연구원.
E-mail : sjlee1110@daum.net
김 형 종(Hyungjong Kim)
2002년 영남대학교 전기전자공학부 졸업.
2005년 한양대학교 전자통신전파공학부 석사 졸업. 2006년∼2009년 ㈜LS산전 주임연구원. 2009년∼현재 서울대학교 전기컴퓨터 공학부 박사과정.
E-mail : hyungjong@gmail.com
남 경 태(Kyung-Tae Nam)
1995년 영남대학교 전자공학과 졸업. 1997년 동 대학원 전자공학과 석사 졸업. 1997년∼1999년 ㈜프로칩스 멀티미디어 주임연구원. 1999년∼2003년 ㈜월텍정보통신 통신연구소 선임연구원. 2003년∼2006년 고등기술연구소 로봇/생산기술센터 선임연구원. 2006년∼현재 한국생산기술연구원 수석연구원.
E-mail : robotnam@kitech.re.kr
국 태 용(Tai-Young Kuc)
1988년 서울대학교 제어계측공학과 졸업.
1993년 포항공과대학 전자전기공학과 박사 졸업. 1995년∼현재 성균관대학교 정보통신공학부 교수.
E-mail : tykuc@yurim.skku.ac.kr
References
Spong M. W. , Hutchinson S. , Vidyasagar M. 2006 Robot Modeling and Control Wiley
Tomei P. 1991 “A simple PD controller for robots with elastic joints” IEEE Trans. Autom. Control 36 (10) 1208 - 1213    DOI : 10.1109/9.90238
Jankovic M. 1995 “Observer based control for elastic joint robots” IEEE Trans. on Robot. Autom. 11 (4)
Bang J. S. , Shim H. , Park S. K. , Seo J. H. 2010 “Robust tracking and vibration suppression for a two-inertia system by combining backstepping approach with disturbance observer,” IEEE Trans. Ind. Electron. 57 (9) 3197 - 3206    DOI : 10.1109/TIE.2009.2038398
Kostarigka A. K. , Doulgeri Z. , Rovithakis G. A. 2013 “Prescribed performance tracking for flexible joint robots with unknown dynamics and variable elasticity,” Automatica 49 1137 - 1147    DOI : 10.1016/j.automatica.2013.01.042
Yun J. N. , Su J. B. 2014 “Design of a disturbance ovserver for a tow link manipulator with flexible joints,” IEEE Trans. Control Syst. Technol. 22 (2) 809 - 815    DOI : 10.1109/TCST.2013.2248733
Moberg S. , Ohr J. , Gunnarsson S. 2008 “A benchmark problem for robust control of a multivariable nonlinear flexible manipulator,” Proc. 17th IFAC World Congr. 1206 - 1211
Moberg S. , Ohr J. , Gunnarsson S. 2009 “A benchmark problem for robust feedback control of a flexible manipulator,” IEEE Trans. Control Syst. Technol. 17 (6) 1398 - 1405    DOI : 10.1109/TCST.2008.2006755
Moberg S. , Wernholt E. , Hanssen S. , Brongardh T. 2014 “Modeling and parameter estimation of robot manipulators using extended flexible joint models,” J. Dyn. Syst. Meas. Control 136
Lotfi N. , Namvar M. 2010 “Global adaptive estimation of joint velocities in robotic manipulators,” IET Control Theory Appl. 4 (12) 2672 - 2681    DOI : 10.1049/iet-cta.2009.0374
Luca A. E. , Schroder D. , Thummel M. 2007 “An acceleration-based state observer for robot manipulators with elastic joints,” IEEE Int. Conf. Robot. Autom. 3817 - 3823
Staufer P. , Gattringer H. 2012 “State estimation on flexible robots using accelerometers and angular rate sensors,” Automatica 22 1043 - 1049
Axelsson P. , Karlssn R. , Norrlof M. 2012 “Bayesian state estimation of a flexible industrial robot,” Automatica 20 1220 - 1228
Chen W. , Tomizuka M. 2012 “Load side state estimation in robot with joint elasticity,” IEEE/ASME Int. Conf. Adv. Intell. Mech. 598 - 603
Chen W. , Tomizuka M. 2014 “Comparative study on state estimation in elastic joints,” Asian J. Control 16 (3) 818 - 829    DOI : 10.1002/asjc.755
Raghavan S. , Hedrick J. K. 1994 “Observer design for a class of nonlinear systems,” Int. J. Control 59 (2) 515 - 528    DOI : 10.1080/00207179408923090
Song B. , Hedrick J. K. 2011 “Nonlinear observer design for lipschitz nonlinear systems,” Proc. Amer. Control Conf. 2578 - 2583
Khalil H. K. 2002 Nonlinear Systems 3rd ed. Prentice Hall
Chen C. T. 1984 Linear System Theory and Design 2nd ed. Oxford University Press
Wonham W. M. 1968 “On a matrix Riccati equation of stochastic control,” SIAM J. Control Optim. 6 (4) 681 - 697    DOI : 10.1137/0306044