Quantification of parameter and model uncertainty for shape memory alloy bending actuators

Abstract

In this article, we employ a Bayesian framework to estimate parameter and model uncertainty for shape memory alloy bending actuators. The Bayesian framework provides parameter densities, instead of ordinary least-squares optimal point estimates. Bayes’ rule relates a posterior parameter density to a prior density and likelihood. However, the posterior density is difficult to calculate directly for high-dimensional parameter spaces. Markov chain Monte Carlo methods overcome this difficulty indirectly by creating a Markov chain whose stationary density is the posterior. In this article, we utilize the Delayed Rejection Adaptive Metropolis algorithm for estimating parameter uncertainty. The shape memory alloy bending actuator is modeled using the homogenized energy framework, a computationally efficient and accurate model for various transductive materials. The model is summarized, and techniques for estimating the heat transfer parameters are presented. An algorithmic approach to quantifying uncertainty is useful for numerous reasons. The anticipated use is to quantify uncertainty for robust control algorithms. Robust control is an area of considerable research for smart materials such as shape memory alloy; however, the source of uncertainty is rarely quantified. The methods employed here would greatly aid in the design of robust controllers.

Keywords

Shape memory uncertainty quantification Markov chain Monte Carlo

Introduction

Shape memory alloys (SMAs) are unique actuators that are gaining increased utilization in prototypes such as robotic catheters (Crews and Buckner, 2012; Veeramani et al., 2008), robotic hands (Laurentis and Mavroidis, 2002), jet chevrons (Hartl et al., 2010a, 2010b), and underwater vessels (Garner et al., 2000). SMA actuators are capable of recovering large strains (approximately 5%) upon heating. However, the design, optimization, and control of SMA devices are complicated by the material’s nonlinear, hysteretic dependence on stress and temperature.

This complex behavior has motivated significant research in the area of robust control, such as sliding mode (or variable structure) control (Ashrafiuon and Jala, 2009; Elahinia and Ashrafiuon, 2002; Elahinia et al., 2005; Esfahani and Elahinia, 2010; Hannen et al., 2011; Shaw et al., 2011; Song et al., 2003; Song and Quinn, 2004; Tai and Ahn, 2011) and $H_{\infty}$ control (Choi et al., 2001; Jayender et al., 2008). The goal of these controllers is to guarantee stable tracking performance in the presence of model uncertainty. However, the authors often fail to ascertain the source of the uncertainty or quantify it, which is the primary purpose for adopting a robust control framework.

In this article, we present a computationally efficient model for a SMA bending actuator that utilizes the homogenized energy model (HEM). Furthermore, we provide a systematic approach to quantifying uncertainty in the model parameters and model output using Markov chain Monte Carlo (MCMC) methods. The Bayesian framework is based on Bayes’ rule, which relates the posterior density of model parameters to a prior density and the likelihood of those parameters (Gelman et al., 2004). However, the posterior density is difficult to calculate for high-dimensional parameter spaces. MCMC methods overcome this difficulty by creating a Markov chain whose stationary density is the posterior density. Here, we use the Delayed Rejection Adaptive Metropolis (DRAM) algorithm (Haario et al., 2006). The advantage of this approach (and all Bayesian methods) is that we obtain densities for model parameters, instead of optimal point estimates as provided in ordinary least-squares (OLS) estimation. It is possible to estimate parameter densities using OLS approaches and the delta method; however, the approach assumes that the distributions are Gaussian. MCMC methods are able to estimate arbitrary densities.

An algorithmic approach to quantifying uncertainty is useful for numerous reasons. As noted, the anticipated use for SMA actuators (in particular the bending actuator under consideration) is to quantify model uncertainty for control algorithms. Robust control requires some measure of model uncertainty in order to guarantee stability or optimality. For example, sliding mode control requires bounds on model uncertainty in order to ensure the attractiveness of the sliding surface (Khalil, 2002). Utilizing an algorithmic approach to parameter uncertainty quantification would greatly aid in the design of these controllers. Additionally, the results can be used to produce bounds on model output for comparison to experimental data and for model prediction with quantified uncertainties.

The remainder of this article is organized as follows: The model of a flexible structure actuated by a single SMA tendon is presented first. The HEM is used, which has been applied to numerous transductive materials (Crews et al., 2012b; Hu et al., 2011, 2012, in press; Smith, 2005; Smith et al., 2003, 2005a, 2005b, 2006). Next, data-driven techniques for estimating model parameters are summarized. The approach to quantify uncertainty is described, including optimization of model parameters and the MCMC method used here. The experimental setup is briefly described. Finally, the results are presented, including parameter densities and correlations and the forward propagation of model uncertainty.

HEM of a single-tendon SMA bending actuator

The system consists of a flexible structure actuated by a single-SMA tendon, as shown in Figure 1(a). The SMA tendon (actuator) is held a fixed distance from the neutral axis of the flexible structure using rapid-prototyped collets. The SMA actuator is prestrained before it is attached to the axially stiff, laterally compliant flexible backbone. Therefore, a moment is created as the SMA tendon contracts due to Joule heating and the structure bends (Figure 1(b)). Similar robotic systems include catheters and smart inhalers (Furst et al., 2010; Furst and Seelecke, 2011). A mesoscopic free energy model of a catheter was introduced in the study by Veeramani et al. (2008), with further studies by Crews (2011) and Crews and Buckner (2012).

System model

As detailed in the studies by Crews (2011) and Crews and Buckner (2012), the bending angle $θ (t)$ can be related to the SMA tendon stress $σ (t)$ by

θ (t) = \frac{a A_{c} L σ (t)}{EI}

(1)

where $a$ is the actuator offset from the neutral axis, $A_{c}$ is the cross-sectional area of the SMA actuator, $L$ is its length, $E$ is the elastic (Young’s) modulus of the flexible structure, and $I$ is its area moment of inertia. The stress–strain relationship for the bending system can be simplified to

σ (t) = \frac{EI}{a^{2} A_{c}} (ε_{P} - ε (t))

(2)

where $ε_{P}$ is the prestrain in the SMA actuator. Equation (2) reduces the planar bending problem to a one-dimensional problem consisting of a SMA actuator in parallel with a linear spring with stiffness

K = \frac{EI}{a^{2} A_{c}}

(3)

The SMA strain $ε (t)$ is modeled using the HEM. The HEM accounts for material inhomogeneities and interaction effects to accurately and efficiently quantify the macroscopic SMA behavior. A complete derivation of the model and data-driven techniques for estimating model parameters are provided in the study by Crews et al. (2012b). In the study by Crews et al. (2012b), a single SMA actuator in constant–temperature conditions was considered. The relevant relationships and equations are summarized here to tailor the HEM to the bending system under consideration.

The SMA strain depends on the stress $σ (t)$ and temperature $T (t)$ . The HEM quantifies the macroscopic strain

ε (t) = \int_{0}^{\infty} \int_{- \infty}^{\infty} ν_{R} (σ_{R}) ν (σ_{I}) \bar{ε} (σ (t) + σ_{I}, T (t); σ_{R}) d σ_{I} d σ_{R}

(4)

by assuming that the relative stress $σ_{R}$ and interaction stress $σ_{I}$ are manifestations of underlying densities $ν_{R} (σ_{R})$ and $ν_{I} (σ_{I})$ . Different forms can be assumed for the densities, subject to certain conditions (Smith, 2005). Here, we assume that the densities

ν_{R} (σ_{R}) = \frac{1}{C_{1}} \sum_{m = 1}^{M_{α}} \sum_{k = 1}^{K_{α}} α_{k, m} ϕ_{k, m} (σ_{R}), C_{1} = \sum_{m = 1}^{M_{α}} \sum_{k = 1}^{K_{α}} α_{k, m}

(5)

and

ν_{I} (σ_{I}) = \frac{1}{C_{2}} \sum_{k = 1}^{K_{β}} β_{k} ψ_{k} (σ_{I}), C_{2} = \sum_{k = 1}^{K_{β}} β_{k}

(6)

are linear combinations of log-normal and normal densities. The kernels $ϕ_{k, m} (σ_{R})$ and $ψ_{k} (σ_{I})$ are given by

ϕ_{k, m} (σ_{R}) = \frac{1}{c_{k} σ_{R} \sqrt{2 π}} \exp (- \ln (σ_{R}) - {μ_{R_{m}}}^{2} / 2 c_{k}^{2})

and

ψ_{k} (σ_{I}) = \frac{1}{b_{k} \sqrt{2 π}} \exp (- \frac{σ_{I}^{2}}{2 b_{k}^{2}})

The density coefficients $α_{k, m}$ and $β_{k}$ and parameters $c_{k}$ , $μ_{R_{m}}$ , and $b_{k}$ can be estimated and optimized using experimental data (Crews et al., 2012b).

The mesoscopic (or local) strain $\bar{ε} (σ (t) + σ_{I}, T (t); σ_{R})$ in equation (4) is given by

\begin{matrix} \bar{ε} (σ (t) + σ_{I}, T (t); σ_{R}) = x_{A} (t) \frac{σ (t)}{E_{A}} \\ + x_{M +} (t) (\frac{σ (t)}{E_{M}} + ε_{T}) + x_{M -} (t) (\frac{σ (t)}{E_{M}} - ε_{T}) \end{matrix}

where $x_{A} (t)$ , $x_{M +} (t)$ , and $x_{M -} (t)$ are the austenite ( $A$ ), martensite plus ( $M +$ ), and martensite minus ( $M -$ ) phase fractions, respectively; $E_{A}$ is the austenitic elastic modulus; $E_{M}$ is the martensitic elastic modulus; and $ε_{T}$ is the maximum recoverable strain.

The evolution of the phase fractions is governed by the coupled differential equations

\begin{matrix} {\overset{\cdot}{x}}_{M +} (t) = - p_{+ A} x_{M +} (t) + p_{A +} x_{A} (t) \\ {\overset{\cdot}{x}}_{M -} (t) = - p_{- A} x_{M -} (t) + p_{A -} x_{A} (t) \end{matrix}

(7)

and the conservation relation

x_{A} (t) = 1 - x_{M +} (t) - x_{M -} (t)

(8)

Substituting equation (8) into equation (7) yields

\begin{matrix} {\overset{\cdot}{x}}_{M +} (t) = - (p_{+ A} + p_{A +}) x_{M +} (t) - p_{A +} x_{M -} (t) + p_{A +} \\ {\overset{\cdot}{x}}_{M -} (t) = - (p_{- A} + p_{A -}) x_{M -} (t) - p_{A -} x_{M +} (t) + p_{A -} \end{matrix}

(9)

The transition rates $p_{α β}$ ( $α, β = A, M +, M -$ ) depend on the SMA stress, temperature, and material parameters such as $σ_{R}$ and are described in greater detail in the study by Crews et al. (2012b). The temperature in the SMA actuator is quantified using a first-order, lumped capacitance model

\overset{\cdot}{T} (t) = - h (T (t) - T_{\infty}) + H ({\overset{\cdot}{x}}_{M +} (t) + {\overset{\cdot}{x}}_{M -} (t)) + γ v^{2} (t)

(10)

where $h$ is the convective heat transfer coefficient between the actuator and ambient air at temperature $T_{\infty}$ , $H$ is the latent heat of phase transformation, $v (t)$ is the applied voltage, and $γ$ is a parameter that accounts for the SMA resistance and specific heat.

The integral (equation (4)) is discretized using four-point Gaussian quadrature on 20 equal intervals, yielding

ε (t) = \sum_{i = 1}^{N_{i}} \sum_{j = 1}^{N_{j}} ν_{R} (σ_{R_{i}}) ν_{I} (σ_{I_{j}}) \bar{ε} (σ (t) + σ_{I_{j}}, T (t); σ_{R_{i}}) v_{i} w_{j}

(11)

where $σ_{R_{i}}$ and $σ_{I_{j}}$ are the quadrature points and $v_{i}$ and $w_{j}$ are the quadrature weights. The summations in equation (11) can be expressed as

ε (t) = V^{T} Γ W

where

\begin{matrix} V^{T} = v_{1} ν_{R} (σ_{R_{1}}), \dots, v_{N_{i}} ν_{R} (σ_{R_{N_{i}}}) \\ W^{T} = w_{1} ν_{I} (σ_{I_{1}}), \dots, w_{N_{j}} ν_{I} (σ_{I_{N_{j}}}) \\ Γ = X_{A} \frac{σ (t)}{E_{A}} + X_{M +} (\frac{σ (t)}{E_{M}} + ε_{T}) + X_{M -} (\frac{σ (t)}{E_{M}} - ε_{T}) \end{matrix}

The $N_{i} \times N_{j}$ matrices $X_{A}$ , $X_{M +}$ , and $X_{M -}$ compile the phase fractions evaluated at the quadrature points

\begin{matrix} X_{α} {(t, σ (t), T (t))}_{ij} = x_{α} (t, σ (t) + σ_{I_{j}}, T (t); σ_{R_{i}}) \\ α = A, M +, M - \end{matrix}

where $x_{α}$ are the solutions to equations (7) and (8).

The ODEs (equation (9)) are discretized and solved using an implicit Euler scheme. For discretized time $t^{k} = k Δ t$ , stress $σ^{k} = σ (t^{k})$ , and temperature $T^{k} = T (t^{k})$ , the implicit Euler discretization yields

\begin{matrix} a_{11}^{k + 1} x_{M +}^{k + 1} + a_{12}^{k + 1} x_{M -}^{k + 1} = x_{M +}^{k} + a_{13}^{k + 1} \\ a_{21}^{k + 1} x_{M +}^{k + 1} + a_{22}^{k + 1} x_{M -}^{k + 1} = x_{M -}^{k} + a_{23}^{k + 1} \end{matrix}

(12)

where

\begin{matrix} a_{11}^{k + 1} = 1 + Δ t (p_{+ A} + p_{A +}) & a_{21}^{k + 1} = p_{A -} Δ t \\ a_{12}^{k + 1} = p_{A +} Δ t & a_{22}^{k + 1} = 1 + Δ t (p_{- A} + p_{A -}) \\ a_{13}^{k + 1} = p_{A +} Δ t & a_{23}^{k + 1} = p_{A -} Δ t \end{matrix}

(13)

Note that in equation (13), the transition rates depend on $σ^{k}$ , $T^{k}$ , $σ_{I}$ , and $σ_{R}$ . The austenitic phase fraction $x_{A}^{k + 1}$ is given by the conservation relation (equation (8)). The implicit Euler discretization on the temperature yields

T^{k + 1} = d_{1}^{k + 1} T^{k} + d_{2}^{k + 1}

where

\begin{matrix} d_{1}^{k + 1} = \frac{1}{1 + h Δ t} \\ d_{2}^{k + 1} = \frac{Δ t (h T_{\infty} + H ({\overset{\cdot}{x}}_{M +}^{k} + {\overset{\cdot}{x}}_{M -}^{k}) + γ {(v^{k})}^{2})}{1 + h Δ t} \end{matrix}

The stress $σ^{k}$ is found by substituting equation (11) into equation (2) and solving for the equilibrium stress, which yields

σ^{k} = \frac{ε_{P} - ε_{T} (V^{T} X_{M +} W - V^{T} X_{M -} W)}{\frac{a^{2} A_{c}}{EI} + V^{T} X_{A} W \frac{1}{E_{A}} + V^{T} X_{M +} W \frac{1}{E_{M}} + V^{T} X_{M -} W \frac{1}{E_{M}}}

Equation (12) can be solved using Cramer’s rule, which gives

\begin{matrix} x_{M +}^{k + 1} = c_{11}^{k + 1} x_{M +}^{k} + c_{12}^{k + 1} x_{M -}^{k} + c_{13}^{k + 1} \\ x_{M -}^{k + 1} = c_{21}^{k + 1} x_{M +}^{k} + c_{22}^{k + 1} x_{M -}^{k} + c_{23}^{k + 1} \end{matrix}

where

\begin{matrix} c_{11}^{k + 1} = \frac{a_{22}^{k + 1}}{det} c_{12}^{k + 1} = - \frac{a_{12}^{k + 1}}{det} c_{13}^{k + 1} = \frac{1}{det} (a_{22}^{k + 1} a_{13}^{k + 1} - a_{12}^{k + 1} a_{23}^{k + 1}) \\ c_{21}^{k + 1} = - \frac{a_{21}^{k + 1}}{det} c_{22}^{k + 1} = \frac{a_{11}^{k + 1}}{det} c_{23}^{k + 1} = \frac{1}{det} (a_{11}^{k + 1} a_{23}^{k + 1} - a_{13}^{k + 1} a_{21}^{k + 1}) \\ det = a_{11}^{k + 1} a_{22}^{k + 1} - a_{12}^{k + 1} a_{21}^{k + 1} \end{matrix}

(14)

The computational efficiency of the model can be increased by storing four-dimensional (4D) arrays $P_{α β}$ of the transition rates evaluated at the quadrature points $σ_{R_{i}}$ and $σ_{I_{j}}$ and uniformly distributed values of the stress $σ^{ℓ}$ and temperature $T^{m}$ , where $σ_{min} \leq σ^{ℓ} \leq σ_{max}$ and $T_{min} \leq T^{m} \leq T_{max}$ . For example, the array $P_{+ A}^{ℓ m}$ would have components

{[P_{+ A}^{ℓ m}]}_{ij} = p_{+ A} (σ^{ℓ} + σ_{I_{j}}, T^{m}; σ_{R_{i}})

During implementation, the indices $ℓ$ and $m$ corresponding to values closest to $σ^{k}$ and $T^{k}$ are found, and the $N_{i} \times N_{j}$ matrices $P_{α β}^{ℓ m}$ are used to calculate the matrices corresponding to equation (14), $C_{11}^{k + 1}, \dots, C_{23}^{k + 1}$ . These matrices can be rapidly calculated using point-wise operations. The matrices of the phase fractions can then be updated in one step using

\begin{matrix} X_{M +}^{k + 1} = C_{11}^{k + 1} . \times X_{M +}^{k} + C_{13}^{k + 1} . \times X_{M -}^{k} + C_{13}^{k + 1} \\ X_{M -}^{k + 1} = C_{21}^{k + 1} . \times X_{M +}^{k} + C_{23}^{k + 1} . \times X_{M -}^{k} + C_{23}^{k + 1} \\ X_{A}^{k + 1} = 1 - X_{M +}^{k + 1} - X_{M -}^{k + 1} \end{matrix}

(15)

where componentwise multiplication ( $. \times$ ) and summation are used. In equation (15), 1 is a $N_{i} \times N_{j}$ matrix with all entries equal to 1. A complete overview of the algorithm is summarized in Algorithm 1.

Data-driven techniques for estimating model parameters

The model contains numerous parameters that can be estimated using experimental data. Accurate initial parameter estimates greatly reduce the computation time required during the optimization step. Furthermore, accurate initial estimates may help the optimization algorithms avoid local minima. The majority of the SMA model parameters are estimated using constant–temperature tensile test data (Crews et al., 2012b). Other parameters such as $a$ , $L$ , $ε_{P}$ , and $EI$ are known or specified during fabrication. Here, we describe a method for estimating the heat transfer model parameters $h$ and $γ$ . The initial estimate of the latent heat $H$ is taken from published values (Huang et al., 2012).

To estimate the heat transfer model parameters, a step voltage is applied to the SMA actuator (Figure 2(a)) until a steady-state bending angle is reached (Figure 2(b)). The voltage is then set to 0 and the decay in the bending angle is measured. To estimate the convection coefficient $h$ , we assume that the decay in bending angle is proportional to the decay in temperature (Figure 2(c)). Note that temperature is not measured. The assumption that the decay rates are the same yields

\frac{T (t_{1}) - T (t_{2})}{T (t_{2}) - T (t_{3})} = \frac{θ_{1} - θ_{2}}{θ_{2} - θ_{3}}

(16)

Since the applied voltage is 0 during this period, the temperature decays according to

T (t) = T_{0} \exp (- ht)

(17)

neglecting the latent heat. Substituting equation (16) into equation (17) and using equal time steps $Δ t$ yield

\frac{\exp (- h Δ t) - \exp (- 2 h Δ t)}{\exp (- 2 h Δ t) - \exp (- 3 h Δ t)} = \frac{θ_{1} - θ_{2}}{θ_{2} - θ_{3}}

Finally, solving for $h$ gives

h = \frac{1}{Δ t} \ln (\frac{θ_{1} - θ_{2}}{θ_{2} - θ_{3}})

(18)

After determining $h$ and the SMA model parameters listed in Table 1 using techniques provided in the study by Crews et al. (2012b), we can estimate $γ$ . As the SMA actuator is heated, the equilibrium stress is the intersection of the unloading curves and the load line of the SMA actuator–flexible structure system, as shown in Figure 3. The load line has a slope equal to the equivalent spring stiffness (equation (3)). Using the steady-state bending angle $θ_{ss}$ , we can calculate the stress using equation (1). The steady-state stress $σ_{ss}$ relates to the temperature by

σ_{ss} = σ_{L} + Δ σ_{T} (T_{ss} - T_{L}) - {\bar{σ}}_{R}

(19)

where $σ_{L}$ , $Δ σ_{T}$ , $T_{L}$ , and ${\bar{σ}}_{R}$ are SMA model parameters described in greater detail in the study by Crews et al. (2012b). Using equation (19), we calculate the approximate temperature $T_{ss}$ in the actuator corresponding to the applied voltage $V_{ss}$ . Finally, assuming steady-state conditions, we have

γ = \frac{- h (T_{ss} - T_{\infty})}{V_{ss}^{2}}

(20)

Table 1.

Model parameters and estimation techniques.

Variable	Description	Estimation technique
$α_{k, m}$	Relative stress density coefficients	Crews et al. (2012b)
$β_{k}$	Interaction stress density coefficients	Crews et al. (2012b)
$E_{A}$	Elastic modulus of austenite	Crews et al. (2012b)
$E_{M}$	Elastic modulus of martensite	Crews et al. (2012b)
$σ_{L}$	Martensite transition stress at temperature $T_{L}$	Crews et al. (2012b)
$T_{L}$	Lower transition temperature	Crews et al. (2012b)
$Δ σ_{T}$	Hysteresis loop’s temperature dependence	Crews et al. (2012b)
$ε_{T}$	Maximum recoverable strain	Crews et al. (2012b)
${\bar{σ}}_{R}$	Mean stress hysteresis loop width	Crews et al. (2012b)
$τ$	Relaxation time	Crews et al. (2012b)
$V$	Layer volume	Crews et al. (2012b)
$h$	Convection coefficient	Equation (18)
$γ$	Heat transfer parameter	Equation (20)
$H$	Latent heat of phase transformation	Huang et al. (2012)
$ε_{P}$	SMA actuator prestrain	Determined during fabrication
$a$	SMA actuator offset from the neutral axis	Determined during fabrication
$L$	Flexible backbone length	Determined during fabrication
$EI$	Backbone elastic modulus and area moment of inertia	Determined during fabrication

SMA: shape memory alloy.

A complete list of the model parameters and estimation techniques is provided in Table 1.

Uncertainty quantification framework

Numerous techniques exist to estimate uncertainty in model parameters. Frequentist approaches involve multiple measurements or samples of experimental data. However, collecting large sets of data can be expensive and time consuming. Alternatively, one could fit statistical moments (such as mean and variance) between model and experimental data (Swiler et al., 2008). Here, we employ a Bayesian framework using MCMC methods (Chib and Greenberg, 1995; Green and Mira, 2001; Haario et al., 2001, 2006). First, we fit the model parameters to experimental data using OLS approaches. Then, we use the OLS results to initialize the DRAM algorithm developed by Haario et al. (2006) and available in the study by Laine (2007).

OLS fit of model parameters

The optimization algorithm minimizes the sum of squared error (SSE)

F (\vec{q}) = \frac{1}{2} \sum_{i = 1}^{N} {(\hat{θ} (t_{i}) - θ (t_{i}; \vec{q}))}^{2}

(21)

between the experimental data $\hat{θ} (t_{i})$ and the model $θ (t_{i}; \vec{q})$ . A complete list of the model parameters $\vec{q}$ and their bounds is provided in Table 2. Note that the parameters are scaled using the relation

\vec{q} = \vec{p} . \times \vec{s}

(22)

where $\vec{p}$ is the vector of physical (unscaled) parameters and $\vec{s}$ is the vector of parameter magnitudes.

Table 2.

Model parameters and associated bounds.

Variable	Lower bound	Upper bound	Units
$α_{k, m}$	0	1	–
$β_{k}$	0	1	–
$E_{A}$	20.0	90.0	GPa
$E_{M}$	20.0	50.0	GPa
$σ_{L}$	200	350	MPa
$Δ σ_{T}$	7.0	12.0	MPa/K
$ε_{T}$	3.5	6.00	%
$τ$	0.001	2.0	s
$V$	1.0 × 10⁻²⁸	1.0 × 10⁻²³	m³
$h$	0.1	6.0	1/s
$γ$	0.5	20.0	K/s V²
$H$	10	200	K
$ε_{P}$	1.0	5.00	%
$a$	0.1	3.0	mm
$L$	60	120	mm
$EI$	0.1	5.0	N cm²

Here, we are including all the SMA model parameters in the optimization routine. Alternatively, one could optimize the SMA model parameters (such as $α_{k, m}$ through $V$ in Table 2) using tensile test data (Crews et al., 2012b) and then only optimize the heat transfer and bending model parameters using the measured bending angle. Including all the parameters will produce a lower SSE but neglects the coupling between parameters. The bounds are chosen to keep the parameters within a reasonable range of the parameters fit to the tensile test data in the study by Crews et al. (2012b) or based on fabrication tolerances.

MCMC methods

Bayes’ rule

π (\vec{q} | \vec{y}) = \frac{p (\vec{y} | \vec{q}) π_{pr} (\vec{q})}{\int_{R^{d}} p (\vec{y} | \vec{q}) π_{pr} (\vec{q}) d \vec{q}}

(23)

relates a posterior density $π (\vec{q} | \vec{y})$ to a prior density $π_{pr} (\vec{q})$ , where $p (\vec{y} | \vec{q})$ is the likelihood of observing $\vec{y}$ given parameters $\vec{q}$ . For relatively few parameters, the posterior density can be calculated directly by numerically integrating the denominator of equation (23). However, in higher-dimensional parameter spaces, this approach is infeasible.

MCMC methods avoid the computationally intensive integral in equation (23) by creating a Markov chain whose stationary density is the posterior density $π (\vec{q} | \vec{y})$ . Various MCMC algorithms exist, including the popular Metropolis–Hastings algorithm (Chib and Greenberg, 1995; Roberts and Rosenthal, 2001; Solonen, 2006) and adaptive algorithms (Haario et al., 2001, 2006; Vihola, 2011). Here, we summarize the DRAM algorithm (Haario et al., 2006), which combines Delayed Rejection (DR) (Green and Mira, 2001) with the Adaptive Metropolis (AM) algorithm (Haario et al., 2001).

The standard Metropolis–Hastings algorithm proposes new parameters ${\vec{q}}^{new}$ from previous parameters ${\vec{q}}^{old}$ using a proposal distribution $f (\cdot | {\vec{q}}^{old})$ . For example, one may use a Gaussian distribution

\begin{matrix} f ({\vec{q}}^{new} | {\vec{q}}^{old}) = \frac{1}{{(2 π)}^{\frac{d}{2}} | V |^{\frac{1}{2}}} \\ \exp (- 0.5 {({\vec{q}}^{new} - {\vec{q}}^{old})}^{T} V^{- 1} ({\vec{q}}^{new} - {\vec{q}}^{old})) \end{matrix}

where $d$ is the dimension of $\vec{q}$ and $V$ is the covariance matrix. The new parameters are accepted with probability

α ({\vec{q}}^{old}, {\vec{q}}^{new}) = min {1, \frac{π ({\vec{q}}^{new} | \vec{y}) f ({\vec{q}}^{old} | {\vec{q}}^{new})}{π ({\vec{q}}^{old} | \vec{y}) f ({\vec{q}}^{new} | {\vec{q}}^{old})}}

(24)

In the special case where $f ({\vec{q}}^{old} | {\vec{q}}^{new}) = f ({\vec{q}}^{new} | {\vec{q}}^{old})$ is symmetric (known as the Metropolis algorithm), equation (24) simplifies to

α ({\vec{q}}^{old}, {\vec{q}}^{new}) = min {1, \frac{π ({\vec{q}}^{new} | \vec{y})}{π ({\vec{q}}^{old} | \vec{y})}}

(25)

Assuming Gaussian model errors, the likelihood is given by

p (\vec{y} | \vec{q}) = \frac{1}{{(2 π σ^{2})}^{1 / 2}} \exp (- 0.5 σ^{- 2} SS (\vec{q}))

(26)

where $SS (\vec{q})$ is the SSE, and $σ$ is the measurement error. Note that the fraction in equation (25) reduces to

\frac{π ({\vec{q}}^{new} | \vec{y})}{π ({\vec{q}}^{old} | \vec{y})} = \exp (- 0.5 σ^{- 2} (SS ({\vec{q}}^{new}) - SS ({\vec{q}}^{old})))

The measurement error $σ$ may be approximated using the relation

σ^{2} = \frac{SS ({\vec{q}}^{*})}{N_{i} - d}

where ${\vec{q}}^{*}$ are the OLS optimal parameters obtained by solving equation (21), $N_{i}$ is the number of measurements, and $d$ is again the dimension of $\vec{q}$ . Alternatively, the measurement error can be sampled using the MCMC algorithm. Using the assumption that the measurement errors are normally distributed and the fact that the likelihood (equation (26)) is in the inverse-gamma family, the conjugate prior is

π_{pr} (σ) \propto {(σ^{2})}^{- (α + 1)} \exp (\frac{β}{σ^{2}})

and the resulting posterior density is

π (σ | \vec{q}, \vec{y}) \propto {(σ^{2})}^{- (α + 1 + n / 2)} \exp (- (β + SS (\vec{q}) / 2) / σ^{2})

Therefore, the measurement error $σ_{k}^{2}$ is sampled from an inverse-gamma distribution

σ_{k}^{2} ~ Γ^{- 1} (α + \frac{N_{i}}{2}, β + \frac{SS ({\vec{q}}^{k - 1})}{2})

or its equivalent representation

σ_{k}^{2} ~ Γ^{- 1} (\frac{n_{0} + N_{i}}{2}, \frac{n_{0} σ_{k - 1}^{2} + SS ({\vec{q}}^{k - 1})}{2})

where the parameter $n_{0}$ is usually small ( $< 1$ ). For a more detailed derivation, see Gelman et al. (2004). Finally, assuming uniform (or noninformative) priors $π_{pr} ({\vec{q}}^{old}) = π_{pr} ({\vec{q}}^{new}) = b - a$ on an interval $[a b]$ , equation (25) further reduces to

\begin{matrix} α ({\vec{q}}^{old}, {\vec{q}}^{new}) \\ = min {1, \exp (- 0.5 σ^{- 2} (SS ({\vec{q}}^{new}) - SS ({\vec{q}}^{old})))} \end{matrix}

(27)

For the SMA bending actuator, the intervals correspond to the bounds listed in Table 2 and used in the OLS estimation.

The Metropolis algorithm uses a Gaussian proposal

{\vec{q}}^{new} = {\vec{q}}^{old} + R^{T} \vec{z}

where $\vec{z}$ is a random vector sampled from the standard normal distribution, and $R$ is the Cholesky decomposition of the covariance matrix $V$ . The adaptation updates the covariance

V_{k} = {\begin{matrix} s_{d} cov ({\vec{q}}^{0}, \dots, {\vec{q}}^{k - 1}) + s_{d} ε I & if k mod k_{a} = 0 \\ V_{k - 1} & else \end{matrix}

(28)

based on the previous parameter chains ${\vec{q}}^{0}, \dots, {\vec{q}}^{k - 1}$ every $k_{a}^{th}$ iteration. In equation (28), $s_{d}$ is a scaling parameter, taken as $s_{d} = 2 . 4^{2} / d$ , where $d$ is the dimension of the parameter space (Gelman et al., 2004). The parameter $ε$ and identity matrix $I$ ensure that $V_{k}$ is not singular; however, the default in the study by Laine (2007) is $ε = 0$ unless this produces a singular matrix. The covariance may be initialized as

V_{0} = σ^{2} {(J^{T} ({\vec{q}}^{*}) J ({\vec{q}}^{*}))}^{- 1}

(29)

where $J ({\vec{q}}^{*})$ is the Jacobian of the model error (equation (21)) evaluated using the OLS optimal parameters ${\vec{q}}^{*}$ . However, here we initialize the proposal covariance

V_{0} = [\begin{matrix} {(0.05 | q_{1}^{*} |)}^{2} & 0 \\ ⋱ \\ 0 & {(0.05 | q_{d}^{*} |)}^{2} \end{matrix}]

(30)

as a diagonal matrix based on the magnitude of the OLS optimal parameters. The DRAM algorithm provided in the study by Laine (2007) utilizes equation (30) as the default initial covariance, and heuristically, we have found that equation (30) works well as the initial proposal covariance.

In the standard Metropolis–Hastings and AM algorithms, if ${\vec{q}}^{new}$ is accepted, then

\begin{matrix} {\vec{q}}^{i} = {\vec{q}}^{new} \\ {\vec{q}}^{old} = {\vec{q}}^{new} \\ S S^{old} = S S^{new} \end{matrix}

otherwise

{\vec{q}}^{i} = {\vec{q}}^{old}

However, with DR, if ${\vec{q}}^{new}$ is rejected, another parameter ${\vec{q}}_{2}^{new}$ is proposed from a distribution $f_{2} (\cdot | {\vec{q}}^{new}, {\vec{q}}^{old})$ , where the second proposal may be based on both ${\vec{q}}^{new}$ and ${\vec{q}}^{old}$ or only one. As detailed in the study by Green and Mira (2001), the second proposal is accepted with probability

α_{2} ({\vec{q}}^{old}, {\vec{q}}^{new}, {\vec{q}}_{2}^{new}) = min {1, \frac{π ({\vec{q}}_{2}^{new} | \vec{y}) f ({\vec{q}}_{2}^{new} | {\vec{q}}^{new}) f_{2} ({\vec{q}}_{2}^{new} | {\vec{q}}^{new}, {\vec{q}}^{old}) (1 - α_{1} ({\vec{q}}_{2}^{new}, {\vec{q}}^{new}))}{π ({\vec{q}}^{old} | \vec{y}) f ({\vec{q}}^{old} | {\vec{q}}^{new}) f_{2} ({\vec{q}}^{old} | {\vec{q}}^{new}, {\vec{q}}_{2}^{new}) (1 - α_{1} ({\vec{q}}^{old}, {\vec{q}}^{new}))}}

(31)

which ensures that Markov chain’s stationary density is the posterior density. The DR can be iterated any number of stages, with the ith stage having proposal $f_{i} (\cdot | {\vec{q}}_{i - 1}^{new}, \dots, {\vec{q}}^{new}, {\vec{q}}^{old})$ and acceptance

\begin{matrix} α_{i} ({\vec{q}}^{old}, {\vec{q}}^{new}, \dots, {\vec{q}}_{i}^{new}) = min {1, \frac{π ({\vec{q}}_{i}^{new} | \vec{y}) f ({\vec{q}}_{i}^{new} | {\vec{q}}_{i - 1}^{new}) f_{2} ({\vec{q}}_{i}^{new} | {\vec{q}}_{i - 1}^{new}, {\vec{q}}_{i - 2}^{new}) \dots f_{i} ({\vec{q}}_{i}^{new} | {\vec{q}}_{i - 1}^{new}, {\vec{q}}_{i - 2}^{new}, \dots, {\vec{q}}^{old})}{π ({\vec{q}}^{old} | \vec{y}) f ({\vec{q}}^{old} | {\vec{q}}^{new}) f_{2} ({\vec{q}}^{old} | {\vec{q}}^{new}, {\vec{q}}_{2}^{new}) \dots f_{i} ({\vec{q}}^{old} | {\vec{q}}^{new}, {\vec{q}}_{2}^{new}, \dots, {\vec{q}}_{i}^{new})} \times \\ \frac{(1 - α_{1} ({\vec{q}}_{i}^{new}, {\vec{q}}_{i - 1}^{new})) (1 - α_{2} ({\vec{q}}_{i}^{new}, {\vec{q}}_{i - 1}^{new}, {\vec{q}}_{i - 2}^{new})) \dots (1 - α_{i - 1} ({\vec{q}}_{i}^{new}, {\vec{q}}_{i - 1}^{new}, {\vec{q}}_{i - 2}^{new}, \dots, {\vec{q}}^{old}))}{(1 - α_{1} ({\vec{q}}^{old}, {\vec{q}}^{new})) (1 - α_{2} ({\vec{q}}^{old}, {\vec{q}}^{new}, {\vec{q}}_{2}^{new})) \dots (1 - α_{i - 1} ({\vec{q}}^{old}, {\vec{q}}^{new}, {\vec{q}}_{2}^{new}, \dots, {\vec{q}}_{i - 1}^{new}))}} \end{matrix}

Here, we are using a second-stage DR. If ${\vec{q}}^{new}$ is rejected, the proposal is scaled by

R_{2} = 0.2 R

and the proposal

{\vec{q}}_{2}^{new} = {\vec{q}}^{old} + R_{2} \vec{z}

(32)

is only based on the last accepted values ${\vec{q}}^{old}$ . Using previous assumptions (such as uniform priors), for the second-stage DR, the acceptance ratio (equation (31)) reduces to

α_{2} ({\vec{q}}^{old}, {\vec{q}}_{2}^{new}) = min {1, \exp (- 0.5 σ^{- 2} (SS ({\vec{q}}_{2}^{new}) - SS ({\vec{q}}^{old}))) r_{2} \frac{1 - α_{1} ({\vec{q}}_{2}^{new}, {\vec{q}}^{new})}{1 - α_{1} ({\vec{q}}^{old}, {\vec{q}}^{new})}}

where

r_{2} = \exp (- 0.5 | | ({\vec{q}}^{new} - {\vec{q}}_{2}^{new}) R^{- 1} {| |}^{2} - 0.5 | | ({\vec{q}}^{new} - {\vec{q}}^{old}) R^{- 1} {| |}^{2})

α_{1} ({\vec{q}}_{2}^{new}, {\vec{q}}^{new}) = min {1, \exp (- 0.5 σ^{- 2} (SS ({\vec{q}}_{2}^{new}) - SS ({\vec{q}}^{new})))}

and $α_{1} ({\vec{q}}^{old}, {\vec{q}}^{new})$ is given by equation (27).

The adaptation mechanism (equation (28)) serves as a global adaptive strategy and updates the proposal function accordingly. The DR step provides local adaptation by scaling the proposal when a step leads to rejection. As detailed in the study by Haario et al. (2006), the two strategies complement each other and may lead to faster convergence than the Metropolis algorithm or either strategy used alone. Adaptation helps correct initial proposals where the covariance is too small. However, if the initial covariance is too large, adaptation may proceed slowly. DR helps by scaling down the proposal with each rejected step. A summary of the DRAM algorithm (Haario et al., 2006) for the SMA bending actuator is presented in Algorithm 2.

Results

A prototype of a SMA-actuated flexible structure, as shown in Figure 4, was constructed to gather experimental data. The system consists of a 0.127 mm diameter FLEXINOL SMA actuator (Dynalloy, Inc., Tustin, CA) and a 0.5 mm diameter superelastic Nitinol beam. The bending angle is measured using a trakSTAR three-dimensional (3D) magnetic tracking system (Ascension Technology Corporation, Burlington, VT). Complete details of the experiment are provided in the study by Hannen et al. (2011).

An input voltage consisting of a sinusoidal function, ramp input, and step input of different magnitudes is used to collect the experimental data, as shown in Figure 5(a). This input quantifies the major loop and various minor loops for the SMA actuator. The resulting measured bending angle is shown in Figure 5(b). The data are first used to fit the model parameters using OLS methods. In this case, we use MATLAB’s lsqnonlin function, which uses either a trust–region–reflective algorithm (Coleman and Li, 1994) or the Levenberg–Marquardt algorithm (Marquardt, 1963). The OLS optimal parameters are then used to initialize the DRAM algorithm. However, out of the 16 model parameters fit using OLS methods, only six are sampled with the DRAM algorithm: $a$ , $ε_{P}$ , $h$ , $γ$ , $EI$ , and $H$ . The heat transfer parameters ( $h$ , $γ$ , and $H$ ) are included because it is anticipated that the heat transfer dynamics will be important in robust controllers. Ultimately, the densities provided for the heat transfer parameters will be utilized to derive controller gains. The additional parameters are included because they have a large influence on the bending mechanics and can vary during fabrication. The remaining SMA model parameters are not sampled by the DRAM algorithm.

The DRAM algorithm uses an adaptation interval $k_{a}$ of 100 iterations. The algorithm is run for 50,000 iterations, at which point the results are analyzed for convergence by viewing the burn-in and mixing of the parameter chains. While burn-in and mixing are largely qualitative (graphical) measures of convergence, quantitative measures can be used (Cowles and Carlin, 1996). If the parameter chains have not converged, additional iterations (50,000 at a time) are used. While quantitative measures of convergence are ideal, especially when MCMC analysis of a model requires months of simulation, 50,000 simulations of the SMA bending actuator are accomplished in approximately 1 day, which allows us to conduct a high number of iterations to ensure convergence.

The DRAM parameter chains are shown in Figure 6. After 50,000 iterations, the sample chains did not appear to be well mixed; therefore, the algorithm was run for an additional 50,000 iterations. As shown in Figure 6, the chains are fully burned-in and well mixed in the final 50,000 iterations. Only the final 50,000 iterations are used to calculate parameter densities and means.

Comparisons between the experimental data and the model output using different parameters are shown in Figure 7. Model fits using the initial parameter estimates are shown in Figure 7(a). A comparison between the OLS fit model and measured bending angle is shown in Figure 7(b). The DRAM mean values are used in Figure 7(c). The initial parameter estimates have a SSE of 448.70, the OLS parameters have a SSE of 4.14, and the DRAM mean values have a SSE of 3.78, which demonstrate an additional advantage of MCMC algorithms. In addition to providing parameter densities, MCMC methods theoretically provide globally optimized values. In practice, one may not know whether the global optimum has been reached. One disadvantage of MCMC methods is the extra computational time. A simulation using the input voltage in Figure 5(a) takes about 0.6 s (approximately 0.5 ms per time step). OLS optimization time is measured on the order of tens of minutes, whereas the MCMC method requires approximately a day for 50,000 iterations.

The initial parameter estimates and parameters estimated only with the OLS method are listed in Table 3. Initial, OLS, and DRAM mean values for the parameters that are also sampled by the DRAM algorithm are provided in Table 4. The majority of the OLS and DRAM parameter values are close to the initial estimates, with the exception of the heat transfer parameters.

Table 3.

Comparison between initial estimates of model parameters and ordinary least-squares estimates.

Variable	Initial estimate	Ordinary least-squares value	Units
$E_{A}$	30.7	20.0	GPa
$E_{M}$	26.0	48.3	GPa
$σ_{L}$	295	350	MPa
$Δ σ_{T}$	9.2	7.6	MPa/K
$ε_{T}$	4.36	6.00	%
$τ$	1.81	1.85	s
$V$	1.0 × 10⁻²⁵	8.35 × 10⁻²⁵	m³
$L$	88	96	mm

Table 4.

Comparison between initial estimates of model parameters, ordinary least-squares estimates, and DRAM values.

Variable	Initial estimate	Ordinary least-squares value	DRAM mean value	Units
$a$	1.0	1.75	1.66	mm
$ε_{P}$	4.0	2.38	2.25	%
$h$	0.581	2.26	2.29	1/s
$γ$	2.50	6.80	7.29	K/s V²
$EI$	1.30	1.34	1.29	N cm²
$H$	48.0	101.4	100.3	K

DRAM: Delayed Rejection Adaptive Metropolis.

Parameter densities are calculated using kernel density estimation (KDE) software (Botev, 2007). The densities are shown in Figure 8. One of the advantages of MCMC methods over other uncertainty quantification approaches such as the delta method is that non-Gaussian densities can be identified and estimated. For example, the densities on $ε_{P}$ (Figure 8(b)) and $EI$ (Figure 8(e)) are not symmetric. Additionally, the heat transfer parameter densities (Figure 8(c), (d) and (f)) appear bimodal, although this may be an artifact of the KDE software. While these densities can be used to propagate model uncertainty, we will use them to develop robust control algorithms (Crews et al., 2012a). For example, the 95% credible intervals may be used to determine parameter bounds for sliding mode controllers. It would be useful comparing the intervals obtained using MCMC to the delta method; however, a well-conditioned sensitivity matrix is required to calculate the covariance (equation (29)), which is not the case here.

In addition to providing a global optimum and parameter densities, the DRAM results can also be used to identify correlation in model parameters, as shown in Figure 9. The results indicate that the heat transfer parameters are highly correlated (Figure 9(f), (m) and (n)), and that $a$ and $ε_{P}$ are correlated (Figure 9(a)). These results may be used to identify correlated parameters for model reduction techniques. Reduction techniques improve the computationally efficiency of the model; however, since the HEM is fast enough for real-time control applications, we retain the physical parameters.

The DRAM chains are resampled to propagate the model uncertainty and calculate credible intervals on the model output for input data not used during model validation. The codes available in the study by Laine (2007) provide functions to resample the chains. The resampled chains produce prediction limits, and the measurement noise is added to predictions to produce observation limits. The 95% observation limits for the SMA actuator are shown in Figure 10. As shown in Figure 10, the majority of the experimental data lies within the bounds. Some data still lies outside of the bounds, especially at small bending angles. The results may be degraded by the fact that the residuals for the SMA bending actuator are not independent and identically distributed (i.i.d). Modeling measurement and model error with an autoregressive model may overcome this limitation and lead to better predictions; however, further investigation is necessary.

The forward propagation of model uncertainty reveals that output uncertainty is largest at steady-state values. For example, the step inputs after 120 s have large uncertainty bounds, as large as 15% of the measured bending angle. Propagation of model uncertainty allows one to determine whether the model accurately quantifies the experimental observations. Additionally, the parameter chains can be sampled to make predictions with quantified uncertainty.

Conclusion

In this article, we employed a model of SMA bending actuators that utilize the homogenized energy framework. We then demonstrated the suitability of the DRAM algorithm for estimating parameter uncertainty. The adaptation in the DRAM algorithm efficiently identifies the ideal proposal for sampling the parameter space. We found the method to be better (converge faster) than the standard Metropolis algorithm. The methods presented here provide parameter densities, instead of single optimal values. Furthermore, the homogenized energy framework has been applied to numerous smart materials, and the Bayesian techniques can easily be adapted to those models and materials.

The MCMC methods provide a number of advantages to traditional OLS estimation of model parameters. The main purpose is to identify model and parameter uncertainties. However, the MCMC algorithms can serve as a global optimizer, whereas OLS approaches are sensitive to local minima. The estimation techniques presented here and in the study by Crews et al. (2012b) also aid in fitting model parameters, as initial points may greatly affect the results. MCMC algorithms can also be used to identify correlation between model parameters, which may be important for model reduction techniques.

As smart materials gain increased utilization in robotic applications, we believe the methods presented here will be instrumental for robust control frameworks. The advantages of robust control are well known; however, often robust control research is presented without any indication of how model and parameter uncertainty are quantified. Often, controller gains are simply chosen and tuned. The Bayesian approach to parameter uncertainty would greatly aid in the development of robust controllers. Future study will focus on the feasibility of the approach in the development of robust controllers, and the parameter densities will be used to derive bounds on model uncertainty.

Footnotes

Acknowledgements

The authors would like to thank Dr Greg Buckner and Jennifer Hannen for providing the experimental data used in this article.

Funding

This research was supported in part by the Air Force Office of Scientific Research through the grant AFOSR FA9550-11-10152.

References

Ashrafiuon

Jala

(2009) Sliding mode control of mechanical systems actuated by shape memory alloy. Journal of Dynamic Systems Measurement and Control: Transactions of the ASME 131: 011010.

Botev

(2007) Kernel density estimator. Available at: http://www.mathworks.com/matlabcentral/fileexchange/14034-kernel-density-estimator

Chib

Greenberg

(1995) Understanding the Metropolis–Hastings algorithm. American Statistician 49: 327–335.

Choi

Han

Kim

. (2001) Force tracking control of a flexible gripper featuring shape memory alloy actuators. Mechatronics 11(6): 677–690.

Coleman

(1994) On the convergence of interior-reflective Newton methods for nonlinear minimization subject to bounds. Mathematical Programming 67: 189–224.

Cowles

Carlin

(1996) Markov chain Monte Carlo convergence diagnostics: a comparative review. Journal of the American Statistical Association 91(434): 883–904.

Crews

(2011) Development of a shape memory alloy actuated robotic catheter for endocardial ablation: modeling, design optimization, and control. PhD Thesis, North Carolina State University, Raleigh, NC.

Crews

Buckner

(2012) Design optimization of a shape memory alloy actuated robotic catheter. Journal of Intelligent Material Systems and Structures 23: 545–562.

Crews

Smith

Hannen

(2012a) Development of robust control algorithms for shape memory alloy bending actuators. In: ASME 2012 conference on smart materials, adaptive structures, and intelligent systems, Stone Mountain, 19–21 September 2012, GA.

10.

Crews

Smith

Pender

. (2012b) Data-driven estimation of the homogenized energy model parameters for shape memory alloys. Journal of Intelligent Material Systems and Structures 23(17): 1897–1920.

11.

Elahinia

Ashrafiuon

(2002) Nonlinear control of a shape memory alloy actuated manipulator. Journal of Vibration and Acoustics: Transactions of the ASME 124(4): 566–575.

12.

Elahinia

Ashrafiuon

Ahmadian

. (2005) A temperature-based controller for a shape memory alloy actuator. Journal of Vibration and Acoustics 127: 285–291.

13.

Esfahani

Elahinia

(2010) Developing an adaptive controller for a shape memory alloy walking assistive device. Journal of Vibration and Control 16(13): 1897–1914.

14.

Furst

Seelecke

(2011) Experimental validation of different methods for controlling a flexible nozzle using embedded SMA wires as both positioning actuator and sensor. Proceedings of SPIE 7978.

15.

Furst

Hangekar

Seelecke

(2010) Practical implementation of resistance feedback measurement for position control of a flexible smart inhaler nozzle. In: ASME 2010 conference on smart materials, adaptive structures, and intelligent systems, 28 September – 1 October 2010 Philadelphia, PA.

16.

Garner

Wilson

Lagoudas

. (2000) Development of a shape memory alloy actuated biomimetic vehicle. Smart Materials and Structures 9: 673–683.

17.

Gelman

Carlin

Stern

. (2004) Bayesian Data Analysis. 2nd ed. Chapman & Hall/CRC, Boca Raton, FL.

18.

Green

Mira

(2001) Delayed rejection in reversible jump Metropolis-Hastings. Biometrika 88: 1035–1053.

19.

Haario

Laine

Mira

. (2006) DRAM: efficient adaptive MCMC. Statistics and Computing 16: 339–354.

20.

Haario

Saksman

Tamminen

(2001) An adaptive Metropolis algorithm. Bernoulli 7: 223–242.

21.

Hannen

Crews

Buckner

(2011) Indirect intelligent sliding mode control of a shape memory alloy actuated flexible beam using hysteretic recurrent neural networks. Smart Materials and Structures. 21: 085015.

22.

Hartl

Lagoudas

Calkins

. (2010a) Use of a Ni60Ti shape memory alloy for active jet engine chevron application: I. Thermomechanical characterization. Smart Materials and Structures 19: 1–14.

23.

Hartl

Lagoudas

Calkins

. (2010b) Use of a Ni60Ti shape memory alloy for active jet engine chevron application: II. Experimentally validated numerical analysis. Smart Materials and Structures 19: 1–14.

24.

Smith

Ernstberger

(2012) Data driven techniques to estimate parameters in a rate-dependent ferromagnetic hysteresis model. Physica B 407(9): 1394–1398.

25.

Smith

Ernstberger

(2012) The homogenized energy model (HEM) for characterizing polarization and strains in hysteretic ferroelectric materials: implementation algorithms and data-driven parameter estimation techniques. Journal of Intelligent Material Systems and Structures. 23(16): 1869–1894.

26.

Smith

Stuebner

. (2011) Statistical parameter estimation for macro fiber composite actuators using the homogenized energy model. Proceedings of SPIE 7978.

27.

Huang

Leary

Ataalat . (2012) Optimisation of Ni-Ti shape memory alloy response time by transient heat transfer analysis. Materials & Design 35: 655–663.

28.

Jayender

Patel

Nikumb

. (2008) Modeling and control of shape memory alloy actuators. IEEE Transactions on Control Systems Technology 16(2): 279–287.

29.

Khalil

(2002) Nonlinear Systems. Prentice Hall. Upper Saddle River, NJ.

30.

Laine

(2007) MCMC toolbox for MATLAB. Available at: http://helios.fmi.fi/~lainema/mcmc/

31.

Laurentis

Mavroidis

(2002) Mechanical design of a shape memory alloy actuated prosthetic hand. Technology and Health Care 10: 91–106.

32.

Marquardt

(1963) An algorithm for least-squares estimation of nonlinear parameters. Journal of the Society for Industrial and Applied Mathematics 11(2): 431–441.

33.

Roberts

Rosenthal

(2001) Optimal scaling for various Metropolis-Hastings algorithms. Statistical Science 16: 351–367.

34.

Shaw

Pan

Huang

(2011) Adaptive fuzzy sliding mode control of X–Y platform driven by shape memory alloy. Journal of Intelligent Material Systems and Structures 22(16): 1803–1810.

35.

Smith

(2005) Smart Material Systems: Model Development. Philadelphia, PA: SIAM.

36.

Smith

Seelecke

Ounaies

. (2003) A free energy model for hysteresis in ferroelectric materials. Journal of Intelligent Material Systems and Structures 14(11): 719–739.

37.

Smith

Hatch

Mukherjee

. (2005a) A homogenized energy model for hysteresis in ferroelectric materials: general density formulation. Journal of Intelligent Material Systems and Structures 16: 713–732.

38.

Smith

Seelecke

Dapino

. (2005b) A unified framework for modeling hysteresis in ferroic materials. Journal of the Mechanics and Physics of Solids 54: 46–85.

39.

Smith

Dapino

Braun

. (2006) A homogenized energy framework for ferromagnetic hysteresis. IEEE Transactions on Magnetics 42(7): 1747–1769.

40.

Solonen

(2006) Monte Carlo methods in parameter estimation of nonlinear models. PhD Thesis, Lappeenranta University of Technology, Finland.

41.

Song

Quinn

(2004) Experimental study of the robust tracking control of a shape memory alloy wire actuator. Journal of Dynamic Systems Measurement and Control: Transactions of the ASME 126: 674–677.

42.

Song

Chaudhry

Batur

(2003) Precision tracking control of shape memory alloy actuators using neural networks and a sliding-mode based robust controller. Smart Materials and Structures 12: 223.

43.

Swiler

Adams

Eldred

(2008) Model calibration under uncertainty: matching distribution information. In: Proceedings of the 12th AIAA/ISSMO. Multidisciplinary Analysis and Optimization Conference, 10–12 September 2008, Victoria, British Columbia, AIAA-2008-5944.

44.

Tai

Ahn

(2011) Adaptive proportional–integral–derivative tuning sliding mode control for a shape memory alloy actuator. Smart Materials and Structures 20: 055010.

45.

Veeramani

Buckner

Owen

. (2008) Modeling the dynamic behavior of a shape memory alloy actuated catheter. Smart Materials and Structures 19(1): 1–14.

46.

Vihola

(2011) On the stability and ergodicity of adaptive scaling Metropolis algorithms. Stochastic Processes and Their Applications 121: 2839–2860.