Online Inverse Optimal Control for Time-Varying Cost Weights

Abstract

Inverse optimal control is a method for recovering the cost function used in an optimal control problem in expert demonstrations. Most studies on inverse optimal control have focused on building the unknown cost function through the linear combination of given features with unknown cost weights, which are generally considered to be constant. However, in many real-world applications, the cost weights may vary over time. In this study, we propose an adaptive online inverse optimal control approach based on a neural-network approximation to address the challenge of recovering time-varying cost weights. We conduct a well-posedness analysis of the problem and suggest a condition for the adaptive goal, under which the weights of the neural network generated to achieve this adaptive goal are unique to the corresponding inverse optimal control problem. Furthermore, we propose an updating law for the weights of the neural network to ensure the stability of the convergence of the solutions. Finally, simulation results for an example linear system are presented to demonstrate the effectiveness of the proposed strategy. The proposed method is applicable to a wide range of problems requiring real-time inverse optimal control calculations.

1. Introduction

The integration of biological principles with robotic technology heralds a new era of innovation, with a significant focus on applying optimal control and optimization methods to analyze animal motion. This approach guides robotic movement development evident in [1], which explores the intricate control systems in mammalian locomotion. Such research underpins the development of robots that emulate the efficiency and adaptability found in nature.

These advancements in understanding animal locomotion through optimal control methods set the stage for the relevance of inverse optimal control (IOC). IOC offers a retrospective analysis of expert movements—human or animal—to infer underlying cost functions optimized in these motions. This methodology is crucial when direct modeling of optimal strategies is complex or unknown.

The use of inverse optimal control (IOC) to identify suitable cost functions from the observable control input and state trajectories of experts is becoming increasingly important. Several successful applications of IOC in estimating the cost weights of multi-features have been reported. For example, the knowledge and expertise of specialists can be categorized and exploited in several fields, including robot control and autonomous driving. The authors of [2], who employed game theory in tailoring robot–human interactions, proposed a method for estimating the human cost function and selecting the robot’s cost function based on the results, leading to the Nash equilibrium in human–robot interactions. The authors of [3] applied IOC to analyze taxi drivers’ route choices. To investigate the cost combination of human motion, the authors of [4] conducted an experiment using IOC techniques to study human motion during the performance of a goal-achieving task using one arm. Additionally, the authors of [5] represented the learning of biological behavior as an inverse linear quadratic regulator (LQR) problem and proposed adaptive methods for modeling and analyzing human reach-to-grasp behavior. Furthermore, the authors of [6] employed an IOC method to segment human movement.

Linear quadratic regulation is a common optimal control method for linear systems. In the 1960s and 1970s, numerous researchers offered solutions to the inverse LQR problem [7,8,9]. Recently, the theory of linear matrix inequality was employed to solve the inverse LQR problem [5,10,11]. Regarding the application of the IOC method for nonlinear systems, several approaches involving methods such as passivity-based condition monitoring [12] or robust design [13] have been reported.

Recent studies in the field of IOC have demonstrated significant advancements. The authors of [14] provided a comprehensive review of the methodologies and applications in inverse optimization, highlighting its growing importance across various domains. The authors of [15] introduced a novel method for sequential calculation in discrete-time systems, enhancing the IOC model’s efficacy under noisy data conditions. The authors of [16] employed a multi-objective IOC approach to explore motor control objectives in human locomotion, which has implications for predictive simulations in rehabilitation technology. Furthermore, the authors of [17] delve into cost uniqueness in quadratic costs and control-affine systems, shedding light on the non-uniqueness cases in IOC. Moreover, a recent thesis [18] introduces a Collage-Based Approach for solving unique inverse optimal control problems, leveraging the Collage method for ODE inverse problems in conjunction with Pontryagin’s Maximum Principle.

Feature-based IOC methods, which involve modeling the cost function as a linear combination of various feature functions with unknown weights, have gained acclaim in recent years [19,20,21,22]. However, it may be difficult to apply these methods to the analysis of complex, long-term behaviors using simple feature functions, e.g., analyzing human jumping [23]. To address this challenge, the authors of [24] proposed a technique for recovering phase-dependent weights that switch at unknown phase-transition points. This method employs a moving window along the observed trajectory to identify the phase-transition points, with the window length determined by a recovery matrix aimed at minimizing the number of observations required for successful cost-weight recovery. Although this method is effective in estimating phase-dependent cost weights, the complex computational requirements limit its use in real-time applications, such as human–robot collaboration tasks. Additionally, in this method, the cost weights in each phase are assumed to be fixed, which may not be generalizable. For example, the human jump motion in [23] was analyzed using time-varying, continuous cost weights.

Overall, the IOC still has several shortcomings that need to be addressed, particularly when applied in approximating complex, multi-phase, continuous cost functions in real time. In this paper, we propose a method for recovering the time-varying cost weights in the IOC problem for linear continuous systems using neural networks. Our approach involves constructing an auxiliary estimation system that closely approximates the behavior of the original system, followed by determining the necessary conditions for tuning the weights of the neurons in the neural network to obtain a unique solution for the IOC problem. We demonstrate that the unique solution corresponds to achieving a zero error between the original system state and the auxiliary estimated system state, as well as zero error between the original costate and the integral of the estimated costate. Based on this analysis, we develop two neural-network frameworks: one for approximating the cost-weight function and the other for addressing the error introduced by the auxiliary estimation system. Additionally, we discuss the necessary requirements for the feature functions to ensure the well-posedness of our online IOC method. Finally, we validate the effectiveness of our method through simulations.

This work makes several significant contributions:

• We provide a solution for the recovery of time-varying cost weights, essential for analyzing real-world animal or human motion.
• Our method operates online, suitable for a broad spectrum of real-time calculation problems. This contrasts with previous online IOC methods that mainly focused on constant cost weights for discrete system control.
• We introduce a neural network and state observer-based framework for online verification and refinement of estimated cost weights. This innovation addresses the critical need for solution uniqueness and robustness against data noise in IOC applications.

2. Problem Formulation

2.1. System Description and Problem Statement

Consider an object’s system dynamics formulated as

$$\dot{x}=Ax+Bu$$

(1)

where $A\in R^{n\times n}$ and $B\in R^{n\times m}$ are two time-invariant matrixes, $x\in R^n$ represents the system states, and $u={[u_1,\dots ,u_m]}^T\in R^m$ denotes the control input of the system [25].

To minimize the following cost function while accounting for dynamics (1), the classic optimal control problem is required to design the optimal control input $u^{*}\left(t\right)$, and generate a sequence of optimal states $x^{*}\left(t\right)$. (Superscript ∗ stands for the optimal condition.)

$$V(x,t)={\int }_t^{t_f}{L}_0(x,u,\tau )d\tau$$

(2)

Here, $L_0$ has the following form:

$$L_0=q^{T}F\left(x\right)+r^{T}G\left(u\right)$$

(3)

where $q={[q_1,q_2,\dots ,q_{n_f}]}^T\in R^{n_f}$ and $r={[r_1,r_2,\dots ,r_m]}^T\in R^m\forall r_i>0$ represent the cost weight vectors, $F\left(x\right)$ is referred to as the general union feature vector with respect to x, and $G\left(u\right)$ indicates the feature vector that is only relevant to the control input u [26]. $n_f$ represents the feature’s number, which is different from the dimension of system states. For simplicity, we assume that $r^{T}G\left(u\right)=u^{T}Ru$ where R is an unknown matrix with $R=\left[\begin{array}{ccc}{r}_1& 0& \cdots \\ ⋮& \ddots & ⋮\\ \cdots & 0& r_m\end{array}\right]$. Additionally, it is assumed that ($A,B$) is controllable, B is a full column rank matrix, and A and B are bounded such that $\left|\right|A\left|\right|\le {\delta }_A\left|\right|B\left|\right|\le {\delta }_B$.

2.2. Maximum Principle in Forward Optimal Control

To minimize the cost function as is the case in (2) with $L_0$ defined in (3), there exists a costate variable vector $\lambda$ that satisfies Pontryagin’s maximum principle as follows:

$$\begin{array}{c}\hfill \dot{\lambda }=-{\overline{F}}_x^{T}q-A^T\lambda \end{array}$$

(4)

$$\begin{array}{c}\hfill Ru+B^T\lambda =0\end{array}$$

(5)

where ${\overline{F}}_x=\frac{\partial F\left(x\right)}{\partial x}$ and $\lambda \in R^n$ denote the costate variables. These two equations are derived from Pontryagin’s Maximum Principle by taking the partial derivatives of the Hamiltonian function defined by $\mathbb{H}(x,u,\lambda )=L_0+{\lambda }^T(Ax+Bu)$, specifically $\dot{\lambda }=-\frac{\partial \mathbb{H}}{\partial x}$ and $\frac{\partial \mathbb{H}}{\partial u}=0$. The initial value of $\lambda$ can be represented as ${\lambda }_0$.

The optimal control input $u^{*}$ of the system expressed by (1) is given as

$$u^{*}=-R^{-1}{B}^T\lambda$$

(6)

where $\lambda$ is unknown. Thus, using this optimal control input, we have

$$\dot{x}=Ax-H\lambda$$

(7)

where H denotes the matrix $H=B{R}^{-1}{B}^T$. Notably, given that B is a full column rank matrix, it is clear that H is invertible. In addition, since B is a bounded constant matrix, there exists a positive scalar ${\delta }_H$ such that H satisfies $\left|\right|H\left|\right|\le {\delta }_H$.

Additionally, the time derivatives of the system dynamics can be formulated as follows:

$$\ddot{x}=A\dot{x}-H\dot{\lambda }$$

(8)

2.3. Analysis of the IOC Problem

We assume that the system states $x[t,t_f]$ and the control input $u[t,t_f]$, which represent the time series of the system states and control inputs from time point t to $t_f$, provide the solution to the optimal minimization of the cost function (2). In addition, we assume that the optimal system states and control input satisfy the boundary conditions $\left|\right|x\left|\right|\le {\delta }_x$ $\left|\right|u\left|\right|\le {\delta }_u$ $\left|\right|\dot{u}\left|\right|\le {\delta }_{\dot{u}}$.

The objective of the IOC problem is to recover the unknown cost weight’s vector $q\left(t\right)$. Furthermore, IOC, for example, may be employed to analyze different behaviors such as the effect of different occasions on the relative importance of certain human motion feature functions. A rigorous analysis of the derived cost weights that can recreate the original data $x[t,t_f],u[t,t_f]$ is required for the aforementioned applications. To begin, we consider two problems:

• What happens when a different feature function is selected?

In previous studies, it was assumed that the cost weight vector q is either a constant value [19] or a step function with multiple phases [24]. These assumptions have been effective in recovering the cost weights used in the analysis of optimal control methods for a robot’s motion control, such as analyzing the motion of a robot controlled by a LQR approach. However, occasionally, it may be inappropriate to assume that the cost weights are constants or step functions when analyzing the complex behaviors of natural objects, such as human motion. In particular, deciding which feature function to adopt when evaluating the motion of natural objects could pose a challenge.

Proposition 1.

Depending on the different selections of feature functions $F\left(x\right)$ for the IOC, the original constant cost weight q may become a time-varying continuous function.

Proof. 

From (8), for the objects’ original feature function, we have

$$\begin{array}{c}\hfill H^{-1}(-\ddot{x}+A\dot{x}+H{A}^T{H}^{-1}Bu)={\overline{F}}_{ox}^T{q}_o\end{array}$$

(9)

where $q_o$ denotes the original time-invariant cost weight vector, and ${\overline{F}}_o\left(x\right)$ denotes the partial derivative with respect to x of the original feature function. When we choose a different feature function $F_n\left(x\right)$, the above equation becomes

$$\begin{array}{c}\hfill H^{-1}(-\ddot{x}+A\dot{x}+H{A}^T{H}^{-1}Bu)={\overline{F}}_{nx}^T{q}_n\end{array}$$

(10)

where ${\overline{F}}_{nx}$ denotes the partial derivative with respect to x of the new selected feature function and $q_n$ is the corresponding cost weights on ${\overline{F}}_{nx}$. Thus, we have

$${\overline{F}}_{ox}^T{q}_o={\overline{F}}_{nx}^T{q}_n\forall t_0\le t\le t_f$$

From this equation, it follows that $q_n$ may be a time-varying function when ${\overline{F}}_{ox}$ and ${\overline{F}}_{nx}$ are not equivalent, and as ${\overline{F}}_{ox}$ and ${\overline{F}}_{nx}$ are continuous functions, we can reasonably conclude that $q_n$ is also a continuous function.    □

Based on this proposition, it is crucial to expand the definition of cost weights to include time-varying values, as this will facilitate a more accurate analysis of the motion of increasingly complex natural objects. Despite the need for time-varying cost weight recovery in many applications, it has received minimal research attention thus far.

• Whether or not the given set $x[t,t_f],u[t,t_f]$ in the IOC problem has a unique solution $\left\{q\right(t),r\}$.

The uniqueness of the solution to the IOC problem when cost weights are constant has been discussed in many studies [15,17,18,22]. In this work, we determine if there is still a unique solution to the IOC problem when q is a time-varying function.

From (10), we can find different continuous functions $q\left(t\right)$ such that the equation is satisfied for different values of R (different values of H). This implies that if q is considered as a time-varying function, the set $\left\{q\right(t),r\}$ will not have a unique solution.

Therefore, when we consider the unique solution of the IOC problem with the time-varying function $q\left(t\right)$, it is necessary to introduce additional conditions to ensure that the IOC problem has a unique solution and that the resulting unique solution is meaningful.

In this study, for simplicity, we assume that $R=I$ [27,28], where I is the identity matrix. In actual optimal control cost functions, when we focus on reducing one of the control inputs $u_i$, the convergence of the i-th system state $x_i$ related to $u_i$ will also be affected. Consequently, the final control result shows that the change in each state of the system is not solely influenced by the chosen cost weights $q\left(t\right)$, but also by $R\left(t\right)$. In the IOC problem, setting $R\left(t\right)=I$ allows the effect of different weights on different control inputs in the original system to be reflected in the current estimate of $q\left(t\right)$. This enables us to view the estimated weights on the system states as representing the relative importance of each state in the system’s dynamic evolution, without considering the impact of the control input on these weights.

Based on our conclusion that q may be time-varying when different feature functions are chosen and on the corresponding conditions under which a unique solution exists, we can define the IOC problem to be solved in this study as follows:

Problem 1.

Online Estimation of Time-Varying Cost weights $q\left(t\right)$

Given: (1) Measured system state x as well as control input u (2) $R=I$

Goal: Online estimate of the time-varying $q\left(t\right)$ utilizing the given x and u.

3. Adaptive Observer-Based Neural Network Approximation of Time-Varying Cost Weights

In this study, we estimate time-varying cost weight functions online using an observer-based adaptive neural network estimation approach, as opposed to earlier studies that required a large number of time series of x and u to recover fixed cost weights offline.

Construction of the Observer

Following the introduction of $\widehat{q}\left(t\right)\in R^n$ denoting the estimation of $q\left(t\right)$, we define the estimation of the associated costate variable $\widehat{\lambda }$ as follows:

$$\dot{\widehat{\lambda }}=-{\overline{F}}_{\widehat{x}}^T\widehat{q}\left(t\right)-A^T\widehat{\lambda }$$

(11)

where ${\overline{F}}_{\widehat{x}}=\frac{\partial F\left(\widehat{x}\right)}{\partial \widehat{x}}$ denotes the partial derivatives of the feature functions that are only relevant to the estimated system states $\widehat{x}$ obtained by inserting $\widehat{\lambda }$ into (7):

$$\dot{\widehat{x}}=A\widehat{x}-H\widehat{\lambda }$$

(12)

where the initial state ${\widehat{x}}_0$ of this system is selected to be ${\widehat{x}}_0=x_0$.

Thus, compared with that of the original system, the error generated by the new estimation system can be expressed as

$$\dot{\tilde{x}}=A\tilde{x}-H\tilde{\lambda }$$

(13)

$$\dot{\tilde{\lambda }}=-{\overline{F}}_x^{T}q\left(t\right)+{\overline{F}}_{\widehat{x}}^T\widehat{q}\left(t\right)-A^T\tilde{\lambda }$$

(14)

where $\tilde{\lambda }=\lambda -\widehat{\lambda }$ and $\tilde{x}=x-\widehat{x}$. Here, the feature function is selected such that its partial derivative with respect to x is bounded and it is assumed that $\left|\right|{\overline{F}}_x^x\left|\right|\le {\delta }_{nx}$, $\left|\right|{\overline{F}}_{\widehat{x}}\left|\right|\le {\delta }_{n\widehat{x}}$ and $\left|\right|{\overline{F}}_x^x\left(x\right)-{\overline{F}}_{\widehat{x}}\left(\widehat{x}\right)\left|\right|\le \zeta \left|\right|\tilde{x}\left|\right|$ where ${\delta }_{nx}$, ${\delta }_{n\widehat{x}}$ and $\zeta$ denote a positive scalar.

Additionally, the time derivatives of (14) can be expressed as

$$\begin{array}{c}\hfill \ddot{\tilde{x}}=A\dot{\tilde{x}}-H\dot{\tilde{\lambda }}\end{array}$$

(15)

Thus, the following equation can be satisfied:

$$\begin{array}{c}\hfill \dot{s}=A_{r}s+T_x\tilde{q}+(T_x-T_{\widehat{x}})\widehat{q}\end{array}$$

(16)

where $s=\left[\begin{array}{c}\dot{\tilde{x}}\\ \tilde{\lambda }\end{array}\right]$, $A_r=\left[\begin{array}{cc}A& H{A}^T\\ 0& -A^T\end{array}\right]$, $T_x=\left[\begin{array}{c}H{\overline{F}}_x\\ -{\overline{F}}_x\end{array}\right]$, $T_{\widehat{x}}=\left[\begin{array}{c}H{\overline{F}}_{\widehat{x}}\\ -{\overline{F}}_{\widehat{x}}\end{array}\right]$. $\tilde{q}$ denotes the error of estimating q. Here, $\left|\right|{\overline{F}}_x^x\left(x\right)-{\overline{F}}_{\widehat{x}}\left(\widehat{x}\right)\left|\right|\le \zeta \left|\right|\tilde{x}\left|\right|$ implies that there exists a positive scalar ${\zeta }^{\prime }$ such that $\left|\right|T_x-T_{\widehat{x}}\left|\right|\le {\zeta }^{\prime }\left|\right|\tilde{x}\left|\right|$ holds. Based on the bound of ${\overline{F}}_x^x\left(x\right),{\overline{F}}_{\widehat{x}}\left(\widehat{x}\right),H$, it follows that there are two positive scalars ${\delta }_{t_x}$ and ${\delta }_{t_{\widehat{x}}}$ such that the following inequalities hold: $\left|\right|T_x\left|\right|\le {\delta }_{t_x}$ and $\left|\right|T_{\widehat{x}}\left|\right|\le {\delta }_{t_{\widehat{x}}}$.

Moreover, from (6) and (7), $\lambda$ can be calculated as follows:

$$\lambda =-H^{-1}Bu$$

(17)

4. Neural Network-Based Approximation of Time Varying Cost Weights

In this section, a neural network-based cost weight approximation algorithm is proposed. To calculate an approximation of the time-varying vector q, we adopt a neural network in which the chosen inputs are $u_I=\left[\begin{array}{c}{x}_0\\ u\end{array}\right]$, where $x_0$ denotes the initial state of the system (1). Based on this, we assume that time-invariant weight matrixes $W\in R^{n_f\times l}$ exist that satisfy the following expression:

$$\begin{array}{c}\hfill q=W^T\varphi \left(u_I\right)+{ϵ}_1\left(u_I\right)\end{array}$$

(18)

where $\varphi \left(u_I\right)$ denotes the activation function and ${ϵ}_1\left(u_I\right)$ denotes the structure approximation error of the neural networks. In addition, the activation function selected enables the activation function as well as its partial derivative to satisfy the following boundary condition: $\left|\right|\varphi \left(u_I\right)\left|\right|\le {\delta }_p$ and $\left|\right|\frac{\partial \varphi \left(u_I\right)}{\partial u_I}\left|\right|\le {\delta }_{pu}$ where ${\delta }_p$ and ${\delta }_{pu}$ represent two positive scalars. Additionally, $\left|\right|{ϵ}_1\left(u_I\right)\left|\right|\le {ϵ}_n$ where ${ϵ}_n$ is a positive scalar.

The estimate of vector q is constructed as follows:

$$\begin{array}{c}\hfill \widehat{q}={\widehat{W}}^T\varphi \left(u_I\right)\end{array}$$

(19)

where $\widehat{W}$ denotes the estimation of W. In this paper, we will combine two estimators ${\widehat{W}}_1$ and ${\widehat{W}}_2$ to estimate W, as shown in Section 4.1. Before presenting the details of the estimators, we first discuss the necessary conditions for the estimation.

Based on the setting of estimator $\widehat{W}$, the error of estimating q can be expressed as

$$\begin{array}{c}\hfill \tilde{q}=q-\widehat{q}={\tilde{W}}^T\varphi \left(u_I\right)+{ϵ}_1\left(u_I\right)\end{array}$$

(20)

where $\tilde{W}=W-\widehat{W}$ denotes the error of estimating W. Substituting $\tilde{q}$ into (16) yields

$$\begin{array}{c}\hfill \dot{s}=A_{r}s+T_x{\tilde{W}}^T\varphi \left(u_I\right)+(T_x-T_{\widehat{x}})\widehat{q}+T_x{ϵ}_1\left(u_I\right)\end{array}$$

(21)

To profoundly comprehend the necessary condition for the convergence of the estimation error $\tilde{W}$, we define uniformly ultimately bounded (UUB) below.

Definition 1.

A time-varying signal $\sigma \left(t\right)$ can be said as UUB if there exists a compact set $S\subset R^n$ so that for all $\sigma \in S$, there exists a bound $\mu \ge 0$ and a time T such that $\left|\right|\sigma \left|\right|\le \mu$ for all $t\ge t_0+T$.

Lemma 1.

If the following conditions are satisfied, $\tilde{W}$ becomes UUB.

• ${\int }_{t_0}^{t_i}sdt,s$ become UUB after a time point $t_1$ ($\left|\right|{\int }_{t_0}^{t_i}sdt\left|\right|\le {\delta }_1$, and $\left|\right|s\left|\right|\le {\delta }_2$)
• The change in $\widehat{W}$ approaches zero
• Matrix C defined below will become a full row rank matrix.

$$\begin{array}{c}\hfill C=\left[\begin{array}{c}{\int }_{t_1+1}^{t_1+2}{T}_x{(I\otimes \varphi \left(u_I\right))}^{T}dt\\ ⋮\\ {\int }_{t_i-1}^{t_i}{T}_x{(I\otimes \varphi \left(u_I\right))}^{T}dt\end{array}\right]\end{array}$$

(22)

where $t_1\le t_i\le t_f$ and any term in C satisfies the persistent excitation (PE) condition defined below.

$$\begin{array}{c}\hfill \left|\right|{\int }_{t_j}^{t_j+1}{T}_x(I\otimes \varphi \left(u_I\right)){dt)}^T\left|\right|\ge {\beta }_j\forall t_1\le t_j\le t_i\end{array}$$

(23)

Here, ${\beta }_j$ is a positive value.

Proof. 

From (21)

$$\begin{array}{cc}\hfill s& =A_r{\int }_{t_0}^{t_i}sdt+{\int }_{t_0}^{t_i}{T}_x{\tilde{W}}^T\varphi \left(u_I\right)dt\hfill \\ \hfill & +{\int }_{t_0}^{t_i}(T_x-T_{\widehat{x}})\widehat{q}dt+{\int }_{t_0}^{t_i}{T}_x{ϵ}_1\left(u_I\right)dt\hfill \end{array}$$

(24)

Since ${\int }_{t_0}^{t_i}sdt\to 0,s\to 0$ reaches a steady state and $A_r$ is a constant, we can obtain the following:

$$\begin{array}{c}\hfill \left|\right|s-A_r{\int }_{t_0}^{t_i}sdt\left|\right|\le {\delta }_{si}\end{array}$$

(25)

where ${\delta }_{si}$ denotes a small positive scalar. Additionally, with both ${ϵ}_1\left(u_I\right)$ and $T_x$ being bounded, this leads to

$$\begin{array}{c}\hfill \left|\right|{\int }_{t_0}^{t_i}{T}_x{ϵ}_1\left(u_I\right)dt\left|\right|\le {\delta }_{Tϵ}\end{array}$$

(26)

where ${\delta }_{Tϵ}$ denotes a small positive scalar. The term ${\int }_{t_0}^{t_i}{T}_x{ϵ}_1\left(u_I\right)dt$ captures the effect of the structural error of the neural network on state s. Since $T_x$ is bounded, when the neural network approximates the cost weight function adequately, the value of ${ϵ}_1\left(u_I\right)$ decreases, which in turn minimizes the overall integral value. In other words, a well-selected neural network structure with a good approximation of the cost weight function will produce a small structure error and, therefore, a small overall integral value ${\int }_{t_0}^{t_i}{T}_x{ϵ}_1\left(u_I\right)dt$.

(24)–(26) leads to

$$\begin{array}{c}\hfill \left|\right|{\int }_{t_0}^{t_i}{T}_x{\tilde{W}}^T\varphi \left(u_I\right)dt+{\int }_{t_0}^{t_i}(T_x-T_{\widehat{x}})\widehat{q}dt\left|\right|\le {\delta }_{si}+{\delta }_{Tϵ}\end{array}$$

(27)

Similarly, we can obtain a similar relation for the duration $[t_0,t_1]$

$$\begin{array}{c}\hfill \left|\right|{\int }_{t_0}^{t_1}{T}_x{\tilde{W}}^T\varphi \left(u_I\right)dt+{\int }_{t_0}^{t_1}(T_x-T_{\widehat{x}})\widehat{q}dt\left|\right|\le {\delta }_{si}+{\delta }_{Tϵ}\end{array}$$

(28)

From (27) and (28), it follows that   

$$\begin{array}{c}\hfill \left|\right|{\int }_{t_1+1}^{t_i}{T}_x{\tilde{W}}^T\varphi \left(u_I\right)dt+{\int }_{t_1+1}^{t_i}(T_x-T_{\widehat{x}})\widehat{q}dt\left|\right|\le 2({\delta }_{si}+{\delta }_{Tϵ})\end{array}$$

(29)

Furthermore, considering ${\int }_{t_0}^{t_i}sdt\to 0$ after $t_1$, the definition of s and $\left|\right|T_x-T_{\widehat{x}}\left|\right|\le {\zeta }^{\prime }\left|\right|\tilde{x}\left|\right|$, this implies that

$$\begin{array}{cc}\hfill \left|\right|{\int }_{t_1+1}^{t_i}(T_x-T_{\widehat{x}})\widehat{q}dt\left|\right|& \le {\int }_{t_1+1}^{t_i}\left|\right|(T_x-T_{\widehat{x}})\left|\right|\left|\right|\widehat{q}\left|\right|dt\hfill \\ & \le {\int }_{t_1+1}^{t_i}{\zeta }^{\prime }{\delta }_{\tilde{x}}{\delta }_{\widehat{q}}dt\equiv {\delta }_{\zeta (t_i-t_1-1)}\hfill \end{array}$$

(30)

where ${\delta }_{\tilde{x}}$ and ${\delta }_{\widehat{q}}$ represent the bounds of $\tilde{x}$ and $\widehat{q}$ respectively. Thus, this leads to the inequality

$$\begin{array}{c}\hfill \left|\right|{\int }_{t_1+1}^{t_i}{T}_x{\tilde{W}}^T\varphi \left(u_I\right)dt\left|\right|\le 2({\delta }_{si}+{\delta }_{Tϵ})+{\delta }_{\zeta (t_i-t_1-1)}\end{array}$$

(31)

In this case, when $\dot{\widehat{W}}$ approaches zero, the following relation emerges:

$$\begin{array}{cc}\hfill & \left|\right|{\int }_{t_1+1}^{t_i}{T}_x{(I\otimes \varphi \left(u_I\right))}^{T}vec\left(\tilde{W}\right)dt\left|\right|\hfill \\ \hfill & =\left|\right|{\int }_{t_1+1}^{t_i}{T}_x{(I\otimes \varphi \left(u_I\right))}^{T}dtvec\left(\tilde{W}\right)\left|\right|\hfill \\ \hfill & \le 2({\delta }_{si}+{\delta }_{Tϵ})+{\delta }_{\zeta (t_i-t_1-1)}\hfill \end{array}$$

(32)

Based on this relation, it follows that

$$\begin{array}{c}\hfill \left|\right|{\int }_{t_1+1}^{t_1+2}{T}_x{(I\otimes \varphi \left(u_I\right))}^{T}dtvec\left(\tilde{W}\right)\left|\right|\\ \hfill \le 2({\delta }_{si}+{\delta }_{Tϵ})+{\delta }_{\zeta \left(1\right)}\end{array}$$

(33)

where ${\delta }_{\zeta \left(1\right)}={\int }_{t_1+1}^{t_1+2}{\zeta }^{\prime }{\delta }_{\tilde{x}}{\delta }_{\widehat{q}}dt=\cdots ={\int }_{t_i-1}^{t_i}{\zeta }^{\prime }{\delta }_{\tilde{x}}{\delta }_{\widehat{q}}dt$.

Thus, it implies that

$$\begin{array}{c}\hfill \left|\right|Cvec\left(\tilde{W}\right)\left|\right|\le (t_i-t_1-1)(2({\delta }_{si}+{\delta }_{Tϵ})+{\delta }_{\zeta \left(1\right)})\end{array}$$

(34)

where C is defined in (22). Due to C being full row rank, this leads to

$$\begin{array}{cc}\hfill \left|\right|vec\left(\tilde{W}\right)\left|\right|& \le \left|\right|C^{+}\left|\right|\left|\right|Cvec\left(\tilde{W}\right)\left|\right|\hfill \\ \hfill & \le \left|\right|C^{+}\left|\right|(t_i-t_1-1)(2({\delta }_{si}+{\delta }_{Tϵ})+{\delta }_{\zeta \left(1\right)})\hfill \end{array}$$

(35)

From (23), we have $\left|\right|C^{+}\left|\right|\le \frac{1}{\sqrt{(t_i-t_1-1){\beta }_j^2}}$

$$\begin{array}{c}\hfill \left|\right|vec\left(\tilde{W}\right)\left|\right|\le \sqrt{\frac{t_i-t_1-1}{{\beta }_j^2}}(2({\delta }_{si}+{\delta }_{Tϵ})+{\delta }_{\zeta \left(1\right)})\end{array}$$

(36)

Thus, $\tilde{W}$ is UUB.

Notably, ${\beta }_j$ evaluates the lower bound of the norm of ${\int }_{t_j}^{t_j+1}{T}_x(I\otimes \varphi \left(u_I\right)){dt)}^T$, it can increase when the data x cause the norm of the integral to deviate significantly from zero. The size of ${\delta }_{\zeta \left(1\right)},{\delta }_{si}$ is related to the minimization of s and ${\int }_{t_0}^{t_i}sdt$, and the size of ${\delta }_{Tϵ}$ is related to the approximation ability of the chosen neural network. The bound of $\tilde{W}$ after $t_1$ can be minimized by the excited x, successfully minimizing s and ${\int }_{t_0}^{t_i}sdt$ while appropriately designing the structure of the neural network.    □

4.1. Construction of the Neural Network

As shown in Lemma 1, the convergence of ${\int }_{t_0}^{t}sd\tau$ is essential in the convergence of $\tilde{W}$ to 0. Therefore, it is necessary to incorporate this consideration in the approximation design.

First, we divide the estimation of the weights of the neural network into two parts:

$$\widehat{W}={\widehat{W}}_1+{\widehat{W}}_2$$

(37)

and

$$\begin{array}{c}\hfill \widehat{q}={\widehat{q}}_1+{\widehat{q}}_2={({\widehat{W}}_1+{\widehat{W}}_2)}^T\varphi \left(u_I\right)\end{array}$$

(38)

where ${\widehat{q}}_1={\widehat{W}}_1^T\varphi \left(u_I\right)$ and ${\widehat{q}}_2={\widehat{W}}_2^T\varphi \left(u_I\right)$.

The necessity for employing two distinct estimators, ${\widehat{W}}_1$ and ${\widehat{W}}_2$, is rooted in their specialized roles in minimizing the tracking error s. This dual-estimator approach ensures that $\widehat{q}\left(t\right)$ closely aligns with the desired trajectory $q\left(t\right)$. While ${\widehat{W}}_1$’s adaptive tuning is primarily aimed at steering s towards zero, its inherent residual errors in its adaptive process necessitate the deployment of ${\widehat{W}}_2$ for error compensation and enhanced accuracy in tracking the ideal cost weight $q\left(t\right)$. To gain a deeper understanding of this system, we will begin by examining the error dynamics, which forms a fundamental basis for the subsequent detailed exploration of the tuning laws for each estimator.

The state equation describing the error dynamics can be obtained as follows:

$$\dot{s}=A_{r}s+T_x{\tilde{q}}_1+(T_x-T_{\widehat{x}}){\widehat{q}}_1-T_{\widehat{x}}{\widehat{q}}_2$$

(39)

where $s=\left[\begin{array}{c}\dot{\tilde{x}}\\ \tilde{\lambda }\end{array}\right]$, $A_r=\left[\begin{array}{cc}A& H{A}^T\\ 0& -A^T\end{array}\right]$, $T_x=\left[\begin{array}{c}H{\overline{F}}_x\\ -{\overline{F}}_x\end{array}\right]$, $T_{\widehat{x}}=\left[\begin{array}{c}H{\overline{F}}_{\widehat{x}}\\ -{\overline{F}}_{\widehat{x}}\end{array}\right]$.

Further, to effectively minimize ${\int }_{t_0}^{t}sd\tau$, we define vector e as follows:

$$e=(T_x-T_{\widehat{x}}){\widehat{q}}_1+Ks+K_p{\int }_{t_0}^{t}sd\tau -T_{\widehat{x}}{\widehat{q}}_2+A_{r}s$$

(40)

where $K=diag\left([k,\cdots ,k]\right)\in {\mathbb{R}}^{2n\times 2n}$ and $K_p=diag\left([k_p,\cdots ,k_p]\right)\in {\mathbb{R}}^{2n\times 2n}$. Parameters k and $k_p$ are two positive scalars, thus, (39) can be written as:

$$\dot{s}=-Ks-K_p{\int }_{t_0}^{t}sd\tau +T_x{\tilde{q}}_1+e$$

(41)

We suppose that an ideal time-invariant weight matrix $W_2\in R^{n_f\times l}$ exists, which guarantees that

$$\begin{array}{cc}\hfill & (T_x-T_{\widehat{x}}){\widehat{q}}_1+Ks+A_{r}s+K_p{\int }_{t_0}^{t}sd\tau \hfill \\ \hfill & =T_{\widehat{x}}{q}^{\prime }=T_{\widehat{x}}(W_2^T\varphi \left(u_I\right)+{ϵ}_2\left(u_I\right))\hfill \end{array}$$

(42)

where $u_I=\left[\begin{array}{c}{x}_0\\ u\end{array}\right]$.

The estimation error of the neural network can be represented as

$$\begin{array}{cc}\hfill {\tilde{q}}_1& \equiv q-{\widehat{q}}_1={\tilde{W}}_1^T\varphi \left(u_I\right)+{ϵ}_1\left(u_I\right)\hfill \\ \hfill {\tilde{q}}_2& \equiv q^{\prime }-{\widehat{q}}_2={\tilde{W}}_2^T\varphi \left(u_I\right)+{ϵ}_2\left(u_I\right)\hfill \end{array}$$

(43)

and e can be represented as

$$e=T_{\widehat{x}}({\tilde{W}}_2^T\varphi \left(u_I\right)+{ϵ}_2\left(u_I\right))$$

(44)

Therefore, (41) becomes

$$\begin{array}{cc}\hfill \dot{s}& =-Ks-K_p{\int }_{t_0}^{t}sd\tau \hfill \\ \hfill & +T_x({\tilde{W}}_1^T\varphi \left(u_I\right)+{ϵ}_1\left(u_I\right))+T_{\widehat{x}}({\tilde{W}}_2^T\varphi \left(u_I\right)+{ϵ}_2\left(u_I\right))\hfill \end{array}$$

(45)

4.2. Tuning Law of the Neural Network for the Estimation of $q\left(t\right)$

An updating law for a neural network that estimates $q\left(t\right)$ can be represented in Theorem 1, based on the error system’s dynamics that were derived in (45).

Theorem 1.

If we choose the updating laws for the neural network weights ${\widehat{W}}_1$ and ${\widehat{W}}_2$ as shown in (46), respectively, where ${\mathsf{\Gamma }}_1$, ${\mathsf{\Gamma }}_2$, and $k_e$ are positive scalar constants, then state s, ${\int }_{t_0}^{t}sd\tau$ and error e will be UUB.

$$\begin{array}{cc}\hfill {\dot{\widehat{W}}}_1& ={\mathsf{\Gamma }}_1\varphi \left(u_I\right)s^T{T}_x\hfill \\ \hfill {\dot{\widehat{W}}}_2& ={\mathsf{\Gamma }}_2\varphi \left(u_I\right){(s+k_{e}e)}^T{T}_{\widehat{x}}\hfill \end{array}$$

(46)

In addition, if there exist positive constants $t_{\delta }$, ${\beta }_1$, ${\beta }_2$, ${\beta }_3$, and ${\beta }_4$ such that the inequalities in (47) are satisfied for all initial times $t_0$, then the signals ${\tilde{W}}_1$ and ${\tilde{W}}_2$ will also be UUB.

$$\begin{array}{cc}\hfill {\beta }_{2}I& \ge {\int }_{t_0}^{t_0+t_{\delta }}{C}_{p1}{\left(t\right)}^T{C}_{p1}\left(t\right)dt\ge {\beta }_{1}I\hfill \\ \hfill {\beta }_{4}I& \ge {\int }_{t_0}^{t_0+t_{\delta }}{C}_{p2}{\left(t\right)}^T{C}_{p2}\left(t\right)dt\ge {\beta }_{3}I\hfill \end{array}$$

(47)

Here, $C_{p1}\left(t\right)=T_x(I\otimes \varphi {\left(u_I\right)}^T)$, $C_{p2}\left(t\right)=T_{\widehat{x}}(I\otimes \varphi {\left(u_I\right)}^T)$

Proof. 

A proof of this theorem can be found in Appendix A.    □

Applying (46) results in s, ${\int }_{t_0}^{t}sd\tau$, and e being UUB, as shown in Theorem 1. Additionally, (46) shows that when s and e decreases, ${\dot{\widehat{W}}}_1$ and ${\dot{\widehat{W}}}_2$ decrease as well, resulting in a decrease in $\dot{\widehat{W}}={\dot{\widehat{W}}}_1+{\dot{\widehat{W}}}_2$. At this point, as stated in Lemma 1, if the condition of matrix C (defined in Lemma 1), being a full row rank matrix, is satisfied, then $\tilde{W}={\tilde{W}}_1+{\tilde{W}}_2$ will also be UUB. Thus, the solution to the IOC problem can be derived by applying (38).

5. Simulations

5.1. Basic Simulation Conditions

To verify the effectiveness of our method, we performed the simulations using a sample linear system controlled by the optimal control method with the original cost weights R selected in two cases.

The sample linear system dynamics can be formulated as follows:

$$\dot{\theta }=A\theta +B\tau$$

(48)

where $\theta ={[{\theta }_1,{\theta }_2]}^T\in R^2$ represents the system states. We select $A=\left[\begin{array}{cc}30& 80\\ 60& 0\end{array}\right]$, $B=\left[\begin{array}{cc}2& 0\\ 0& 4\end{array}\right]$ and $\tau \in R^2$ denoting the control input.

The cost function selected in these simulations is formulated as

$$\begin{array}{c}\hfill V_r=\frac{1}{2}{\int }_0^{t_f}({\theta }^{T}Q\left(t\right)\theta +{\tau }^{T}R\tau )dt\end{array}$$

(49)

when all the elements of $\theta$ satisfying $|{\theta }_i|\le {\theta }_{r_l}$ and $Q\left(t\right)=\left[\begin{array}{cc}{q}_1& 0\\ 0& q_2\end{array}\right]$ is the continuous time-varying cost weights on system states $\theta$. $R=\left[\begin{array}{cc}{r}_1& 0\\ 0& r_2\end{array}\right]$ represents the cost weights on the control inputs.

Moreover, in our simulations, we select 0 as the initial value of all the elements of both ${\widehat{W}}_1$ and ${\widehat{W}}_2$. Actuation function $\varphi \left(u_I\right)$ was selected as $\varphi \left(u_I\right)={[{\varphi }_1\left(u_I\right),\cdots ,{\varphi }_i\left(u_I\right),\cdots ,{\varphi }_l\left(u_I\right)]}^T$ with ${\varphi }_i\left(u_I\right)$ designed as

$${\varphi }_i\left(u_I\right)=exp\left(\frac{-{(u_I-{\psi }_i)}^T(u_I-{\psi }_i)}{\nu }\right)$$

(50)

where $\nu$ denotes a positive scalar and ${\psi }_i$ denotes the center of the respective activation function. We initialized the activation function centers on a four-dimensional grid to match the dimension of $u_i$, ensuring a uniform distribution across the input space and enhancing network adaptability.

The overall implementation for recovering the time-varying cost weights is shown in Algorithm 1.

Algorithm 1 Online implementation

Input:  $\{x_i,u_i\}$ Output:  $\widehat{q}\left(t\right)$ Initialization: 1: Initialize $\widehat{\lambda }$, $\widehat{x}$, ${\widehat{W}}_1$, ${\widehat{W}}_1$, ${\widehat{W}}_2$, ${\mathsf{\Gamma }}_1$, ${\mathsf{\Gamma }}_2$ and $R=I$. LOOP Process 2: for $i=0$ to K do 3:    Calculate $\lambda$ using $\lambda =-H^{-1}Bu$. 4:    Calculate $\dot{\widehat{x}}$ and $\dot{\widehat{\lambda }}$ using (11) and (12). 5:    Calculate $\dot{\tilde{x}}$ and $\dot{\tilde{\lambda }}$ using (14) and (13). 6:    Calculate $s=\left[\begin{array}{c}\dot{\tilde{x}}\tilde{\lambda }\end{array}\right]$. 7:    Calculate e following (40). 8:    Calculate $\varphi \left(u_I\right)$ and update ${\widehat{W}}_1$, ${\widehat{W}}_2$ using (46). 9:    Calculate $\widehat{q}\left(t\right)$ using (38). 10: end for 11: return$\widehat{q}\left(t\right)$

Two cases are considered in the simulation:

• In the first case, we apply the optimal control of the sample system with cost weights $\theta$ as the signal ($q_1\left(t\right)=1+cos\left(t\right)$ and $q_2\left(t\right)=2+sin\left(t\right)$). The proposed IOC method is employed online to estimate the cost weights, with the simultaneous online recovery of the original system trajectory. Parameters ${\mathsf{\Gamma }}_1$ and ${\mathsf{\Gamma }}_2$ in the updating law are set to ${\mathsf{\Gamma }}_1=1$ and ${\mathsf{\Gamma }}_2=1$, respectively. Parameters k and $k_p$ are set to $k=50$ and $k_p=625$, respectively. The initial values of ${\widehat{W}}_1$ and ${\widehat{W}}_2$ are set to matrixes with all elements equal to zero. The original $r_1$ and $r_2$ are set to $r_1=1$ and $r_2=1$, respectively. The simulation also uses 49 nodes in the neural network.
• In the second case, we perform the simulation of our IOC method, but with the original $r_1$ and $r_2$ set to $r_1=3$ and $r_2=4$, respectively. All other simulation settings are the same as in the first case.

Similar to the simulation sections in previous works ([6,24]), we use the control input from the simulation, which ignores the measurement issues with the control input and measurement errors that may occur in real-world applications. This allows us to purely evaluate the performance of our method in solving the IOC problem. In actual applications, the control input can be calculated by substituting the measured $\dot{\theta }$ into (48), as described in [24].

5.2. Results

The simulation results are shown in the figures below.

In Figure 1, the blue solid line represents the original variation in the cost weights whereas the gray solid line represents the estimated cost weights. After a brief period of oscillation at the initial time, our method accurately recovers the original cost weights when $R=I$. Notably, similar to the case in other adaptive control methods and adaptive neural network based control methods, the initial oscillation is a result of the adaptive initialization of the weights in (46) due to the large initial errors in ${\tilde{W}}_1$ and ${\tilde{W}}_2$.

Figure 2 demonstrates the impact of selecting $R=I$ on the estimation results when the original R value is arbitrary. The solid blue line represents the original time-varying cost weights, whereas the dotted gray line represents the final estimated values. Although the estimated values differ from the original values, the general trend of the changes is preserved. In addition, the gray line represents the mutual weights in the dynamics of the system state, whereas the original weights among the control inputs are reflected in the current estimate of $q\left(t\right)$. From the figure, we can observe that the bottom lines in blue and gray colors represent the value of the original and estimated $q_2$. Evidently, the blue line for $q_2$ is larger than that for $q_1$ from 4.8 s to 5 s. Additionally, in the original settings, $r_2$ is 4, which confers greater importance to the decrease in $u_2$ compared with the case when $r_1=3$, leading to the weakening of the convergence of the ${\theta }_2$ term associated with $u_2$. In our estimates, the value of the dashed line for the estimated $q_2$, which also considers the impact from original setting of R is not greater than the value of estimated $q_1$ between 4.8 s and 5 s. This indicates that the convergence of ${\theta }_2$ is weakened by considering the impact from the cost weights on control input. Our dashed line more accurately reflects the actual situation compared to the blue line.

In Figure 3, Figure 4 and Figure 5, we show the results of error e, states s and ${\int }_{t_0}^{t}sd\tau$ in two cases. The blue lines show the results of the first case, whereas the gray dotted lines show the results of the second case. From the figures, we can observe that all the values effectively decrease to a low range during the simulation, and most importantly, in the second case, the different selections of R do not affect the convergence of these values. This demonstrates the effectiveness of our method and highlights that even with different values of R, the recovered cost weights are still feasible solutions to the IOC problem, as they can be utilized to regenerate a similar system trajectory and control inputs (${\int }_{t_0}^{t}sd\tau =\left[\begin{array}{c}\tilde{x}\\ {\int }_{t_0}^t\tilde{\lambda }d\tau \end{array}\right]\to 0$).

6. Discussion

6.1. Robustness of the Proposed Method to Noisy Data

In (46), ${\mathsf{\Gamma }}_1$ and ${\mathsf{\Gamma }}_2$ decrease the error by regulating the updating speed of the estimated values. Adjusting these two terms may successfully reduce the impact of data noise to a certain degree. Their roles are similar to that of a low-pass filter’s time constant. For example, in the setting of the first case, when noise exists, $x\sim \mathcal{N}(0,{10}^{-1})$ and $u\sim \mathcal{N}(0,{10}^{-4})$, the simulation results show that different sets of ${\mathsf{\Gamma }}_1$ and ${\mathsf{\Gamma }}_2$ (e.g., ${\mathsf{\Gamma }}_1=10$, ${\mathsf{\Gamma }}_2=10$; ${\mathsf{\Gamma }}_1=1$, ${\mathsf{\Gamma }}_2=1$) can significantly influence the noise reduction performance.

As shown in Figure 6, while relatively small values of ${\mathsf{\Gamma }}_1$ and ${\mathsf{\Gamma }}_2$ may result in a low convergence rate, they effectively reduce the impact of data noise. Our method demonstrates robustness against noise by allowing for the adjustment of parameters ${\mathsf{\Gamma }}_1$ and ${\mathsf{\Gamma }}_2$.

6.2. Calculation Complexity and Real-Time Calculation

The proposed algorithm has a low computational complexity, as it only involves the calculation of dot products between matrixes and vectors as well as the summation of vectors. Additionally, it does not require any iterative or optimization calculations. This makes it an efficient solution for real-time calculations. In fact, our simulation shows that a single iteration of the algorithm using case 1 settings takes only approximately 0.23 ms in Matlab 2016b to complete the SIOC’s calculation, which is fast enough to meet real-time calculation requirements.

6.3. Advantages of Using $R=I$

The simulation results suggest that one of the key advantages of setting R as a constant I is that it effectively consolidates the impact of cost weights on state convergence, which would have been influenced by different settings of R, into the estimated value of $q\left(t\right)$. This allows for a comprehensive evaluation of the system state convergence, as it only depends on $q\left(t\right)$, without needing to account for additional considerations. Furthermore, by maintaining a consistent value of $R=I$, it is possible to standardize the analysis of the same motion across multiple agents, which is crucial for various applications.

7. Conclusions

In this paper, we proposed a neural network based method for recovering the time-varying cost weights in the IOC problem for linear continuous systems. Our approach involved constructing an auxiliary estimation system that closely approximates the behavior of the original system, followed by determining the necessary conditions for tuning the weights of the neurons in the neural network to obtain a unique solution for the IOC problem. We discussed the necessary requirements for the previous settings to ensure the well-posedness of our online IOC method. We showed that the unique solution corresponds to achieving a nearly zero error between the original system state and the auxiliary estimated system state, as well as nearly zero error between the original costate and the integral of the estimated costate. Based on this analysis, we developed two neural network frameworks: one for approximating the cost weight function and the other for addressing the error introduced by the auxiliary estimation system and terms. Finally, we validated the effectiveness of our method through simulations, highlighting its ability to recover time-varying cost weights and its robustness against different original choices of R. Overall, our method represents a significant advancement in the field of online IOC, and it is applicable to a wide range of problems requiring real-time IOC calculations.