Distributed Newton method for time-varying convex optimization with backward Euler prediction

1.   Introduction

Time-varying distributed optimization (TVDO) has gained increasing attention due to its practical advantages over time-invariant distributed optimization (TIDO) in dynamic environments where the objective function or constraints may evolve over time. TVDO has been applied in various fields, including power systems ^[1], traffic system ^[2], and robotics ^[3], with further applications in energy management ^[4]. In TVDO, the optimal solution changes over time, making traditional algorithms designed for TIDO, which aim to reach a static optimizer, unsuitable for direct application. To address TVDO, several discrete-time algorithms (DTAs) have been developed ^[5,6,7,8]. For instance, prediction-correction methods are employed in ^[6,7] to solve TVDO with a specific sampling period, where the tracking error is linked to the size of the sampling period. A comprehensive review of DTAs and related work can be found in the survey ^[9]. However, due to factors like sampling period, step size, or errors in local optimization, continuous-time algorithms (CTAs) generally face challenges in asymptotically tracking the time-varying optimal trajectory.

Given the extensive use of digital computing and sensing units ^[10,11], we focus on a discrete-time framework. Specifically, we sample the objective functions $F(x; t)$ at discrete time instances $t_k$, where $k = 0, 1, 2, \dots$ and the sampling period is defined as $h \triangleq t_{k+1}-t_k$. Instead of addressing the time-varying problem directly, we solve a series of time-invariant problems. If the sampling period $h$ is made sufficiently small, the resulting solution trajectory $F(t_k)$ can be obtained with high accuracy. However, in most practical applications, solving these problems for each time sample is impractical, as the computation time required to find each optimizer often surpasses the rate at which the solution trajectory evolves, unless $F(t)$ remains nearly stationary.

The batch method is one traditional approach for addressing the problem sequence, where the objective function is sampled at specified intervals, solving each resultant static problem within given periods. However, this method often fails to align with real-time processing demands, particularly when computational resources are limited and sampling intervals are short. This limitation becomes more pronounced with larger problem sizes, preventing convergence within the required timeframe ^[12]. Consequently, attention has shifted towards online optimization techniques. These methods update the optimization problem continuously throughout the algorithm's iterations, allowing for the extraction of suboptimal solutions at any stage, regardless of convergence status ^[1,13]. Over time, these solutions increasingly approximate the optimal solution. In pursuit of improving this approach, various strategies have been developed. For example, a method based solely on correction operations was introduced, achieving an asymptotic error bound on the order of $O(h)$ ^[14]. Another strategy involves a prediction-correction algorithm tailored for time-varying parameter optimization, though it requires an initial optimal solution, limiting its practical application ^[15]. Nonetheless, these methods have facilitated some theoretical advances in reducing the computational load, especially in convex optimization scenarios. Further, the interior point method has been applied to solve constrained convex optimization problems characterized by time-varying elements, utilizing a log-barrier penalty function and a dynamic system comprising predictive and corrective components based on Newton's method ^[16]. For unconstrained time-varying optimization, a suite of algorithms deploying prediction-correction strategies has been proposed. These include Gradient Trajectory Tracking (GTT) and Approximate Gradient Tracking (AGT), which achieve an $O(h^2)$ asymptotic error range. Moreover, Newton Trajectory Tracking (NTT) and Approximate Newton Tracking (ANT) are shown to offer superior error bounds ^[10]. Building on this framework, a decentralized prediction-correction algorithm for time-varying network optimization challenges has been developed, employing novel matrix splitting techniques and approximating matrix inverses using the Taylor series ^[6,17].

Prediction-correction algorithms ^[18], utilizing nonstationary optimization techniques ^[15,19,20], have been developed to iteratively solve convex programs that change over time. These algorithms function by predicting, at time $t_k$, the optimal solution for the next time step $t_{k+1}$, using an approximation of how the objective function $F$ varies during this interval. The initial prediction is then refined through gradient or Newton descent methods. However, these algorithms are tailored for centralized systems. In our work, we address time-varying convex programs in decentralized environments, where nodes can communicate only with their immediate neighbors. Therefore, the prediction-correction methods proposed in ^[18] are not directly applicable in this decentralized context.

In the discussed algorithms that incorporate a prediction step, while achieving tighter error bounds, these methods necessitate the computation or approximation of the inverse of the Hessian matrix. In certain scenarios, this process also involves calculating mixed partial derivatives. It is established that computing the inverse of a matrix is computationally intensive, with the complexity typically represented as $O(p^3)$, where $p$ is the matrix's dimension. This complexity escalates further for objective functions that are more intricate due to the requirement of mixed partial derivative calculations. In this paper, we outline our major contributions as follows:

1) We introduce a backward Euler prediction step tailored for the unconstrained time-varying optimization problem. This approach eliminates the need for calculating matrix inverses and computing mixed partial derivatives, resulting in reduced computational complexity. The corrective aspect of the algorithm is derived from a Newton step.

2) We establish the convergence of our proposed algorithm and detail its convergence rate. The asymptotic error associated with our algorithm depends on the sampling period $h$, ranging from a worst-case scenario of $O(h^2)$ to an optimal scenario of $O(h^4)$.

The structure of this paper is as follows: In Section Ⅱ, we provide the mathematical foundations essential for the development of the major results. In Section Ⅲ, we detail the algorithms incorporating the backward Euler prediction step. The performance of these algorithms is examined in Section Ⅳ. A numerical example demonstrating the practical application of our theoretical findings is presented in Section Ⅴ. We conclude in Section VI.

2.   Preliminaries

We focus on a connected, undirected graph $\mathcal{G} = (V, E)$, where the vertex set $V$ consists of $n$ nodes, and the edge set $E$ consists of $m$ edges. Distributed optimization algorithms are employed to address the problem of minimizing a global smooth, strongly convex cost function across a set of nodes, where the objective function is expressed as a sum of local functions. We intend to solve the following problem raised in ^[6]

where $\mathbf{x} \in \mathbb{R}^{np}$ represents the stacked vector of decision variables for all nodes, and $\mathbf{x}^i \in \mathbb{R}^p$ denotes the decision variable for node $i$, and $\mathbf{x} = (\mathbf{x}^{1^T}; \dots; \mathbf{x}^{n^T})^T$. The function \(f(\mathbf{x}; t) \colon = \sum\limits_{i \in V} f^i(\mathbf{x}^i; t) \) is the sum of the locally available functions at each node, while \(g(\mathbf{x}; t) \colon = \sum\limits_{i \in V} g^{i, i}(\mathbf{x}^i; t) + \sum\limits_{(i, j) \in E} g^{i, j}(\mathbf{x}^i, \mathbf{x}^j; t) \) incorporates terms induced by the network structure $\mathcal{G}$, necessitating coordination and information exchange across the network. Before introducing distributed protocols to solve the problem in (2.1), we present an example to clarify the problem setting.

For example, in a wireless sensor network, each sensor node must adhere to channel capacity and interference constraints, which can be formulated as a resource allocation problem. The time-varying utility function for sensor $i$ is defined as ${f^i}: \mathbb{R}^p \times \mathbb{R}_+ \to \mathbb{R}$, with the decision variable ${\mathbf{x}^i} \in \mathbb{R}^p$ representing the resources allocated to node $i$ in a network $\mathcal{G}$ consisting of $n$ sensors. This leads to the formulation of the time-varying resource allocation problem.

where the matrix $\mathbf{A} \in \mathbb{R}^{lp \times np}$ is the augmented graph edge incidence matrix. This matrix is composed of $l \times n$ square blocks, each of dimension $p$. For an edge $e = (j, k)$, where $j < k$, that connects node $j$ to node $k$, the block $(e, j)$ in $\mathbf{A}$ is given by $[\mathbf{A}]_{ej} = \mathbf{I}_p$, and the block $[\mathbf{A}]_{ek} = -\mathbf{I}_p$, where $\mathbf{I}_p$ denotes the identity matrix of dimension $p$. All other blocks in $\mathbf{A}$ are zeros. The time-varying vectors $\mathbf{b}(t) \in \mathbb{R}^{lp}$ are determined by channel capacity and rate transmission constraints.

In our work, we consider the nodes to be consumer devices aiming to solve decentralized approximations of the problem defined in (2.2). Specifically, we explore the approximate augmented Lagrangian relaxation of (2.2), and solve

where the parameter $\beta > 0$ acts as a regularization term, encouraging consistency among all nodes. In (2.3), the matrix $\mathbf{A} \in \mathbb{R}^{lp \times np}$ is the graph edge incidence matrix. It is straightforward to observe that the first term in (2.3) is identical to the first term in (2.1). Moreover, given the definition of the augmented graph edge incidence matrix $\mathbf{A}$, we can simplify the second term in (2.3) to ${\left\| {(\mathbf{x}_i - \mathbf{x}_j) - \mathbf{b}_i(t)} \right\|^2}$, which corresponds to the functions ${g^{i, j}}(\mathbf{x}_i, \mathbf{x}_j; t)$ in (2.1).

Notation: Let $\|\cdot\|$ represent the Euclidean norm. The gradient of the function $F(\mathbf{x}; t)$ with respect to $\mathbf{x}$ is denoted by $\nabla_\mathbf{x}F(\mathbf{x}; t)$. The partial derivatives of $\nabla_\mathbf{x}F(\mathbf{x}; t)$ with respect to $\mathbf{x}$ and $t$ are indicated by $\nabla_{\mathbf{xx}}F(\mathbf{x}; t)$ and $\nabla_{t\mathbf{x}}F(\mathbf{x}; t)$, respectively. The third-order derivative of $F(\mathbf{x}; t)$ with respect to $\mathbf{x}$ is denoted as $(\nabla_{\mathbf{xxx}}F(\mathbf{x}; t)$. The time derivative of the Hessian of $F(\mathbf{x}; t)$ with respect to $t$ is expressed as $\nabla_{t\mathbf{xx}}F(\mathbf{x}; t) = \nabla_{\mathbf{x}t\mathbf{x}}F(\mathbf{x}; t)$. Finally, the second derivative of $\nabla_\mathbf{x} F(\mathbf{x}; t)$ with respect to $t$ is represented by $\nabla_{tt\mathbf{x}}F(\mathbf{x}; t)$.

3.   Algorithm

3.1.   The dynamical system

If $\mathbf{x}^*$ represents the optimal solution for the objective function $F(\mathbf{x}; t)$, then the gradient of $F(\mathbf{x}; t)$ with respect to $\mathbf{x}$ must satisfy

As a consequence, the time derivative of this expression must also be zero, leading to

where $\dot{\mathbf{x}}^*(t)$ represents the time derivative of $\mathbf{x}^*(t)$. If the Hessian of $F(\mathbf{x}; t)$ is invertible, we can deduce from the above equation that

We can apply either the gradient descent method or the Newton method to optimize $F(\mathbf{x}; t)$, leading to the formulation of a continuous-time dynamical system as follows:

where $\gamma$ represents a control gain with the constraint $0 < \gamma \le 1$. When $\gamma = \alpha \nabla_{\mathbf{xx}} F(\mathbf{x}; t)$, where $\alpha > 0$ is a constant, this corresponds to the gradient term in Eq (8). Conversely, when $\gamma = 1$, Eq (8) represents the Newton term. The trajectory $\mathbf{x}(t)$ generated by Eq (8) will approach a vicinity of $\mathbf{x}^*(t)$, but due to the time-varying nature of the solution, it does not converge precisely to $\mathbf{x}^*(t)$ ^[21].

If the optimal solution $\mathbf{x}^*(t)$ was known for some $t_0 \geq 0$, the system (3.3) could be used to track the evolution of $\mathbf{x}^*(t)$, since (3.3) guarantees the optimality condition $\nabla_{\mathbf{x}}F(\mathbf{x}^*; t) = 0$ for all $t \geq t_0$. If $\mathbf{x}^*(t)$ is not accessible at any time, we can combine the dynamics (3.3) and (3.4) and devise the following dynamical system:

where $\mathcal{F}(\mathbf{x}; t)$ is defined as

The dynamics described in (3.6) consist of two components: A prediction component and a correction component. The prediction component, given by $-\nabla_{\mathbf{xx}}F(\mathbf{x}; t) \nabla_{t\mathbf{x}}F(\mathbf{x}; t)$, aims to forecast the change in the optimal solution (refer to (3.3)). The correction component, $-\gamma \nabla_{\mathbf{xx}}F^{-1}(\mathbf{x}; t)\nabla_\mathbf{x}F(\mathbf{x}; t)$, serves to steer $\mathbf{x}(t)$ toward the optimum.

3.2.   Discretization calculation

For any time interval $[t_0, t_0 + T]$, where $0 < T < +\infty$, we can partition the interval into $N$ sub-intervals such that $t_{k+1} = t_k + h$, for $k = 0, 1, 2, \ldots, N-1$, where $h$ is the discretization step size. Let $\mathbf{x}(t_k)$ represent the solution of $F(\mathbf{x}; t)$ at time $t_k$. For simplicity, we may denote $\mathbf{x}(t_k)$ as $\mathbf{x}_k$ throughout this paper.

A numerical method for approximating (3.5) is considered a one-step method if, for all $k \geq 0$, $\mathbf{x}_{k+1}$ depends solely on $\mathbf{x}_k$. Two examples of such one-step methods are presented in ^[22].

(1) The forward Euler calculation method

(2) The backward Euler calculation method

The primary distinction between the forward and backward methods is found in their treatment of the first-order approximation term for $F(\mathbf{x}; t)$. It follows from (3.8) that

yielding

Substituting (3.10) into (3.7) gives

According to (3.11), we formulate the prediction step as follows:

The prediction step can be computed only depend on the state information at time $k$ and $k-1$. Compared with the prediction step in ^[6,17], the prediction step in our algorithm can reduce the computation complexity significantly.

To address the time-varying optimization problem in (2.1), we sample the continuous time-varying objective function $F(\mathbf{x}; t)$ at discrete times $t_k$, where $k = 0, 1, 2, \ldots$. This transformation turns the continuous time-varying optimization problem into a sequence of time-invariant optimization problems.

In network, the prediction step can be written as (3.12). For node $i$, the predicted variable $\mathbf{x}^i_{k+1|k}$ is just upon the information of the node itself which given by

the $\mathbf{x}^i_{k+1|k}$ is computed through local operations, but for the correction step in our algorithm, the Newton term involving the calculate of the inverse of the matrix, which requires the global communication. To develop a decentralized protocol to correct the prediction variable $\mathbf{x}_{k+1|k}$ ^[18], here we apply the decentralized Newton method proposed by ^[6]. The decentralized Newton method is formulated as follow

where $\mathbf{c}_{{k+1}, (K)} \in \mathbb{R}^{np}$ is called correction direction, its definition is

For $\mathbf{D}^{-\frac{1}{2}}_{k+1|k}\sum\limits_{\tau = 0}^{K} \left(\mathbf{D}^{-\frac{1}{2}}_{k+1|k}\mathbf{B}_{k+1|k}\mathbf{D}^{-\frac{1}{2}}_{k+1|k}\right)^\tau \mathbf{D}^{-\frac{1}{2}}_{k+1|k}$, it is the K-th order approximate of $\nabla_{\mathbf{xx}}F(\mathbf{x}_{k+1|k}; t_{k+1})^{-1}$, which is derived from truncations of the series ^[6,21]

where matrices $\mathbf{D}_{k+1|k}$ and $\mathbf{B}_{k+1|k}$ are defined as

where the $diag[\nabla_{\mathbf{xx}}g(\mathbf{x}_{k+1|k}; t_{k+1})]$ denotes the block diagonal matrix, which contains the diagonal blocks of the matrix $\nabla_{\mathbf{xx}} g(\mathbf{x}_{k+1|k}; t_{k+1})$. The matrix

is computed at node $i$. The second term in (3.19) links the decisions of node $i$ with those of its neighboring nodes $j \in \mathcal{N}^i$. The structure of the matrix $\mathbf{B}_{k+1|k}$ is determined by the graph topology, with null diagonal blocks $\mathbf{B}^{i, i}_{k+1|k}$ and non-zero off-diagonal blocks $\mathbf{B}^{i, j}_{k+1|k}$ given by $-\nabla_{\mathbf{x}^i\mathbf{x}^j}g^{i, j}(\mathbf{x}^i, \mathbf{x}^j; t)$ whenever nodes $i$ and $j$ are connected. The calculation of $\mathbf{B}^{ij}_{k+1|k}$ is as follows:

Lemma 1. According to ^[21], as per the definition of the correction direction in (3.16) the sequence of correction directions satisfies

define $\mathbf{c}^i_{k+1, (K)}$ and $\nabla_{\mathbf{x}}F^{i}(\mathbf{x}_{k+1|k}; t_{k+1})$ as the $i$-th component of the vector $\mathbf{c}_{k+1, (K)}$ and $\nabla_{\mathbf{x}}F(\mathbf{x}_{k+1|k}; t_{k+1})$, the recursion of (3.21) can be decomposed into local components

The gradient component in (3.22)

is also computed at node $i$, Lemma 1 asserts that the component $\mathbf{c}^i_{k+1, (K)}$ can indeed be calculated using local operations, which means that the iterative application of (3.15) can be executed in a distributed fashion. Consequently, node $i$ computes $\mathbf{x}^i_{k+1}$ by implementing the following local descent:

We call DeNSP as the decentralized Newton method with a backward Euler prediction that uses descent method [cf., (3.24)] and the prediction step [cf., (3.14)]. The algorithm is summarized as Algorithm 1.

4.   Convergence analysis

In this section, we examine the convergence of the algorithms introduced in Section Ⅲ. Our convergence analysis demonstrates that the discrepancy between the optimal solution $\mathbf{x}^*(t_k)$ and the computed solution $\mathbf{x}_k$ is ultimately upper bounded. Specifically, the error is predominantly influenced by the sampling time $h$. Throughout this paper, we adopt the following assumption regarding the objective function:

Assumption 1. Each node's own function $f^i$ is twice differentiable, and the eigenvalues of the local Hessians matrix $\nabla_{\mathbf{x}^i\mathbf{x}^i}f^i(\mathbf{x}^i; t)$ are bounded with positive constants $[m_1, M_1]$, i.e.,

Assumption 2. Functions $g^{i, i}(\mathbf{x}^i; t)$ and $g^{i, j}(\mathbf{x}^i, \mathbf{x}^j; t)$ are twice differentiable and the eigenvalues of the aggregate function Hessian are bounded $[0, L_1]$, the constant $L_1 < \infty$,

Assumption 3. The function $F(\mathbf{x}; t)$ regarding all the derivatives of $\mathbf{x}\in\mathbb{R}^{np}$ and $t\geq0$ are bounded,

From the abounds of Hessians $\nabla_{\mathbf{xx}}f(\mathbf{x}; t)$ and $\nabla_{\mathbf{xx}}g(\mathbf{x}; t)$ in Assumptions 1 and 2 respectively, the Hessian of global cost $\nabla_{\mathbf{xx}}F(\mathbf{x}; t)$ uniformly satisfies

These conditions ensure that the problem formulated in (2.1) is strongly convex, and they also ensure the invertibility of the Hessian matrix. To refine the contribution of our paper, we impose bounds on the higher-order derivatives of $F$, as stated in Assumption 3. A similar assumption has been utilized in previous works ^[6,10,18].

We turn to analyzing the DeNSP algorithm. As per the definition of $\mathbf{D}_{k+1|k}$ and $\mathbf{B}_{k+1|k}$ in (3.18), the matrix $\mathbf{D}_{k+1|k}^{-1/2}\mathbf{B}_{k+1|k}\mathbf{D}^{-1/2}_{k+1|k}$ is positive semidefinite, and its eigenvalues are upper-bounded by a constant $\rho < 1$, as shown in Proposition 2 of ^[17].

where $\rho = (1+\frac{2m_1}{L_1})^{-1} < 1$.

In the following theorem, we establish that the sequence generated by the DeNSP algorithm asymptotically converges to a neighborhood of the optimal trajectory whose radius depends on the sampling period $h$.

Proposition 1. The norm of the difference between prediction $\mathbf{x}_{k+1|k}$ and the optimal solution $\mathbf{x}^*_{k+1}$ is upper bounded by

where $\delta : = 1+h\Big(\frac{D_1(D_0+\gamma D_4)}{m_1^2}+\frac{D_2+\gamma (L_1+M_1)}{m_1} \Big)$, and $\Delta = \frac{D_1(D_0+\gamma D_4)+(D_0+\gamma (L_1+M_1))}{m_1^3}+\frac{(D_0+\gamma D_4)(2D_2+\gamma m_1)}{m_1^2}+\frac{\gamma D_0+D_3}{m_1}$.

Proof. See Supplementary A. □

Theorem 1. Under Assumptions 1–3, fixing $K$ as the level of Hessian inverse approximation for the correction step, there exist bounds $\bar{K}$, $\bar{h}$, and $\bar{R}$, such that if the sampling period $h$ is selected to satisfy $h \leq \bar{h}$, $K$ is chosen to be $K \geq \bar{K}$, and the initial optimality gap fulfills $\| \mathbf{x}_0 - \mathbf{x}^*(t_0)\| \leq \bar{R}$, then $\mathbf{x}_{k+1}$ converges to the optimal trajectory $\mathbf{x}^*(t_{k+1})$ within a bounded error.

specially, when approximation level $K$ is chosen sufficiently large and $\gamma = 1$, then the

Proof. See Appendix B. □

Theorem 1 indicates that the error bound described in (4.7) for the DeNSP algorithm is primarily influenced by the sampling period $h$. In the worst-case scenario, the asymptotic error floor is expected to be on the order of $O(h^2)$. However, according to the result $\lim_{k \to \infty} \sup \left\| \mathbf{x}_{k+1} - \mathbf{x}^*({t_{k+1}}) \right\| = O(h^4)$, it is evident that under certain conditions, the error bound can be further improved. Specifically, if the approximation level $K$ is chosen to be sufficiently large, the DeNSP algorithm has the potential to achieve a much tighter tracking performance, reducing the asymptotic error to the order of $O(h^4)$. This suggests that by carefully selecting the parameter $K$, the algorithm can provide significantly better accuracy in tracking the optimal solution over time. Meanwhile, reducing the sampling interval provides greater assistance in improving convergence accuracy. Furthermore, since the DeNSP algorithm omits the calculation of the Hessian matrix inverse in the prediction step, the computational complexity of the prediction step is $O (p)$, while other algorithms such as DPC-N and DAPC-N require the calculation of the Hessian inverse in the prediction step, ignoring the complexity of gradient calculation, resulting in a computational complexity of $O (p ^ 3)$. In the correction step, the computational complexity of several algorithms is consistent.

5.   Simulation

In this section, we provide a numerical example to demonstrate the effectiveness of the DeNSP algorithm. We consider a resource allocation problem in a network of interconnected devices, as discussed in Section Ⅱ.

We solve the problem described by (2.2) using an approximate augmented Lagrangian method, as outlined in (2.3). The local objective functions associated with sensor $i$ are defined as follows:

In our simulation, we consider the case where decisions are the variables ${\mathbf{x}^i} \in {\mathbb{R}^p}$, $p = 10$, and set the scalar parameters as $w = 0.2\pi$, $k = 0.1$, ${b^{i, l}}\sim\mathcal{U}_{[-2, 2]}^1$ and ${d}(t) = \cos (\omega t)$. We consider there are $n = 50$ nodes in a wireless network which the nodes are randomly distributed in the square $[-1, 1]^2$ and they can only communicate with each other if they are closer than a range of $r = 2.5\sqrt 2 /\sqrt n$. Then the nodes generate a network with $l$ links. We set the $\boldsymbol{b} = 0$ in (2.3) and the tuning parameter $\beta = \sqrt{20}$. We compare the DeNSP algorithm with DPC-N algorithm and DAPC-N algorithm mentioned in ^[6]. The constant step size of DPC-N and DAPC-N algorithm is set to $\gamma = 1$.

From convergence result (4.7), it can be seen that the convergence accuracy is affected by the number of iterations and sampling intervals in the correction step. In order to compare the advantages with other algorithms, we set a fixed and commonly used sampling interval $h = 0.1$. Further observe the impact of the number of iterations in the calibration process on convergence accuracy. In Figure 1, we run different algorithms and show the convergence performance between different algorithms. When the approximation level or the number of communication rounds $K$ is set to 3, the tracking performance of these algorithms is comparable, showing little to no significant difference in their effectiveness. However, as $K$ is increased to 5, the DeNSP algorithm demonstrates a clear advantage, achieving superior tracking accuracy and responsiveness. This suggests that the DeNSP algorithm benefits more from additional communication rounds, which enhances its ability to track the time-varying optimal solution more precisely, making it particularly effective in scenarios where higher precision is required. Further more, the DPC-N algorithm and the DAPC-N algorithm need to calculate of the inverse of the Hessian matrix, with the complexity typically represented as $O(p^3)$, where $p$ is the matrix's dimension, and calculate mixed partial derivatives in the prediction step by communicating with neighbor nodes. In contrast, our algorithm can simply complete the prediction step using the previous information which can reduce computation time greatly.

Although the prediction correction algorithm proposed in this article simplifies the prediction steps, the correction steps still require a large amount of computation on Hessian matrix. Furthermore, for some optimization problems involving integer decision variables, the algorithm proposed in this paper is difficult to solve and requires further research.

6.   Conclusions

In this paper, we propose the backward Euler prediction step for the problem of the unconstrained distributed optimization. Though the theoretical analysis, the convergence accuracy of the proposed algorithm can reach $O(h^2) \sim O(h^4)$. And compared with the DPC-N algorithm and DAPC-N algorithm, our algorithm need not compute the inverse of the Hessian of the cost function, which can save computing time. Finally, we verify the theoretical results via a numerical example. In future research, we will further simplify the calculation of the Hessian matrix in the calibration step.

Author contributions

Zhuo Sun: Conceptualization, Methodology, Writing-review & editing; Huaiming Zhu: Formal analysis, Investigation, Software, Writing-original draft; Haotian Xu: Investigation, Writing-review & editing. All authors have read and approved the final version of the manuscript for publication.

Conflict of interest

The authors declare no conflict of interest.

Supported Foundation

This work is supported in part by the National Natural Science Foundation of China ($61304179, 71431001, 71831002, 72211540399$); the Humanity and Social Science Youth Foundation of the Ministry of Education ($19YJC630151$); the International Association of Maritime Universities ($20200205_AMC$); the Natural Science Foundation of Liaoning Province ($2020-HYLH-32$); the Dalian Science and Technology Innovation Fund ($2020JJ26GX023$)