Interior-point method for NLP

Author: Apoorv Lal (CHEME 6800, Fall 2021)

Introduction

Interior point (IP) methods are used to solve all kind of problems from linear to non-linear, from convex to non convex. The first known IP method is Frisch's (1955) logarithmic barrier method that was later extensively studied by Fiacco and McCormick . However, these methods mainly emerged in the late 1970's and 1980s. They were developed as algorithms to solve linear programming(LP) problems with a better worst case complexity than the simplex method . In 1984, Karmarkar used IP methods to develop 'A new polynomial-Time method for Linear Programming' . The main advantage of this method was better practical performance than the earlier algorithms like the ellipsoid method by Khachiyan . After Karmarkar's publication more intensive research in this field led to development of a new category of IP methods known as 'primal-dual' methods. Currently, most powerful type of primal-dual codes are Ipopt, KNITRO and LOQO . In recent years the majority of research in the field of IPs is for nonlinear optimisation problems, especially for second order (referred as SOCO) and semidefinite optimisation ( referred as SDO). SDO has wide range of applications in combinatorial optimization, control, structural engineering and electrical engineering .

Theory & Methodology

Nonlinear Constrained Optimisation

General form for any nonlinear constrained optimization problem can be depicted as shown in Equation set 1:

$minf(x)$

$s.t.c(x)=0$

$x\geqslant 0$

• The objective function and constraints can be any linear, quadratic or any general nonlinear functions.
• Typically, applied for convex optimization problems
• Convert maximisation problems to minimisation problems by taking negative of the objective function
• Convert any general inequalities to simple inequalities using slack variables.

Slack variables

Consider a general equation of the form as shown in Equation set 2:

$minf(x)$

$s.t.g(x)\geqslant b$

$h(x)=0$

$x\geqslant 0$

Now we need to convert this form to the form as shown in Equation set 1. For this the inequality '$g(x)\geqslant b$ ' can be written $g(x)-b\geq 0$ which can be further written as:

$g(x)-b-s=0$ , along with $s\geq 0$

Thus the set of equation '$g(x)-b-s=0$' and '$h(x)=0$' become '$c(x)=0$' as shown in Equation set 1.

Barrier function

Once the given problem is converted to the standard form (as shown in Equation set 1), we will now try to remove the inequality constraint $x\geqslant 0$ using the barrier term. We will introduce natural log barrier term for inequality constraints as shown in Equation set 3:

$\min f(x)-\mu \textstyle \sum _{i}\displaystyle ln(x_{i})$

$s.t.c(x)=0$

The point $\displaystyle ln(x_{i})$ is not defined for $x_{i}\leqslant 0$ . Our goal is to solve this barrier problem for a given $'\mu \textstyle '$, then decrease this $'\mu \textstyle '$ and resolve the problem. If $'\mu \textstyle '$ is decreased at the right rate and the approximate solutions are good enough then this method will converge to an optimal solution of the non linear optimization problem. We want to only search in the interior feasible space.

Example for barrier function:

Following example shows how barrier term can be incorporated into the objective function:

$min(x-4)^{2}-\mu \textstyle \displaystyle ln(x)$

For different values of $'\mu \textstyle '$ , the value of objective function changes and thus as the value of $'\mu \textstyle '$ is decreased from '10' to '1' (as shown in Figure 1) , we approach towards the optimal solution i.e. $x=4$.

Figure 1: Variation of given objective function for different values of 'μ'

After we know the basics described above, the following methodology adopted from can be used to solve the non-linear optimization problem.

KKT Condition for Barrier Problem

KKT conditions for the barrier problem of the form:

$\min f(x)-\mu \textstyle \sum _{i=1}^{n}\displaystyle ln(x_{i})$ $\Rightarrow$ $\nabla f(x)+\nabla c(x)\lambda -\mu \textstyle \sum _{i=1}^{n}\displaystyle {\tfrac {1}{x_{i}}}$

$s.t.c(x)=0$ ${\displaystyle }$ $c(x)=0$

We can define now define $z_{i}={\tfrac {\mu }{x_{i}}}$ and solve the modified version of the KKT conditions:

$\nabla f(x)+\nabla c(x)\lambda -z=0$

$c(x)=0$

$XZe-\mu e=0$ .

The matrix $'e'$ is simply a column vector of ones i.e. $e={\begin{bmatrix}1\\1\\1\\1\end{bmatrix}}$ to satisfy the rules of matrix operations.

KKT solution with Newton-Raphson method

$\nabla f(x)+\nabla c(x)\lambda -z=0$

$c(x)=0$

$XZe-\mu e=0$

$\Rightarrow$ ${\begin{bmatrix}W_{k}&\nabla c(x_{k})&-I\\\nabla c(x_{k})^{T}&0&0\\Z_{k}&0&X_{k}\end{bmatrix}}$${\begin{pmatrix}d_{k}^{x}\\d_{k}^{\lambda }\\d_{k}^{z}\end{pmatrix}}=$$-{\begin{pmatrix}\nabla f(x_{k})+\nabla c(x_{k})\lambda _{k}-z_{k}\\c(x_{k})\\X_{k}Z_{k}e-\mu _{j}e\end{pmatrix}}$

Here, $d_{k}^{z}$ , $d_{k}^{x}$ and $d_{k}^{\lambda }$ are the search directions for the $x,z,\lambda$ and $W_{k}$ is the second gradient of the Lagrangian.

Thus, $W_{k}={\nabla _{xx}}^{2}L(x_{k},\lambda _{k},z_{k})={\nabla _{xx}}^{2}(f(x_{k})+c(x_{k})^{T}\lambda _{k}-z_{k})$ and $Z_{k}={\begin{bmatrix}z_{1}&0&0\\0&\centerdot \centerdot \centerdot &0\\0&0&z_{n}\end{bmatrix}}$ , $X_{k}={\begin{bmatrix}x_{1}&0&0\\0&\centerdot \centerdot \centerdot &0\\0&0&x_{n}\end{bmatrix}}$

The above system can be arranged into symmetric linear system as follows:

${\begin{bmatrix}W_{k}+\Sigma _{k}&\nabla c(x_{k})\\\nabla c(x_{k})^{T}&0\end{bmatrix}}{\begin{pmatrix}d_{k}^{x}\\d_{k}^{\lambda }\end{pmatrix}}$ $=-{\begin{pmatrix}\nabla f(x_{k})+\nabla c(x_{k})\lambda _{k}\\c(x_{k})\end{pmatrix}}$ where $\Sigma _{k}=X_{k}^{-1}Z_{k}$

We can then solve for $d_{k}^{z}$ after finding the linear solution to $d_{k}^{x}$ and $d_{k}^{\lambda }$ with the following form of the explicit solution: $d_{k}^{z}=\mu _{k}X_{k}^{-1}e-z_{k}-\Sigma _{k}d_{k}^{x}$

Step size $(\alpha )$

After we have found out the search directions we need to find an appropriate value of the step size. The aim here is to find a value of step size to minimise the value of the objective functions and minimise the constraint violations. Following two methods are generally used to determine the value of $(\alpha )$:

• Decreasing the merit function: Merit function is defined as the sum of the objective function and the absolute values of the constraint violations multiplied by a scalar quantity i.e. $merit=f(x)+v\sum \left\vert c(x)\right\vert$ so that it determines if it's a better step or not.
• Filter methods use a combination of objective function $f(x)$ and constraint violations $\left\vert c(x)\right\vert$ to determine whether to accept the step size $(\alpha )$ or not.

Once the step size has been determined we can find out the values of $x_{k+1},z_{k+1},\lambda _{k+1}$

$x_{k+1}=x_{k}+\alpha d_{k}^{x}$

$z_{k+1}=z_{k}+\alpha d_{k}^{z}$

$\lambda _{k+1}=\lambda _{k}+\alpha d_{k}^{\lambda }$

Convergence Criteria

Convergence criteria can be defined when KKT conditions are satisfied within a tolerance:

$\max \left\vert \nabla f(x)+\nabla c(x)\lambda -z\right\vert \leqslant \epsilon _{tol}$

$\max \left\vert c(x)\right\vert \leq \epsilon _{tol}$

$\max \left\vert XZe-\mu e\right\vert \leq \epsilon _{tol}$

Pseudocode

The above described algorithm can be depicted in the form of the following pseudocode as illustrated by

Barrier method vs Primal-dual method

In order to understand the concept of central path and how the barrier method is different from primal-dual method, the following methodology adopted from can be used:

Consider the simple case of an LP with the following form:

$minc^{T}x$

$s.t.Ax\leqslant b$

Numerical Example

Consider the following non-linear optimisation problem:

$\min x_{2}(5+x_{1})$

$s.t.x_{1}x_{2}\geq 5$

$x_{1}^{2}+{x_{2}}^{2}\leq 20$

This problem needs to be converted into standard form first:

$\min x_{2}(5+x_{1})$

$s.t.x_{1}x_{2}-5-x_{3}=0\Rightarrow c_{1}(x)$

$20-x_{1}^{2}-{x_{2}}^{2}-x4=0\Rightarrow c_{2}(x)$

$x_{3}\geq 0$

$x_{4}\geq 0$

The initial guess values can be selected as $x^{o}=$ ${\begin{bmatrix}x_{1}\\x_{2}\\x_{3}\\x_{4}\end{bmatrix}}={\begin{bmatrix}0\\0\\0\\0\end{bmatrix}}$

We will also choose the initial guess values of $\lambda _{1},\lambda _{2},z_{1},z_{2}$ as $'0'$.

$c(x^{o})=$ ${\begin{bmatrix}-5\\20\end{bmatrix}}$

$\nabla _{x}c={\begin{bmatrix}x_{2}&-2x_{1}\\x_{1}&-2x_{2}\\-1&0\\0&-1\end{bmatrix}}={\begin{bmatrix}0&0\\0&0\\-1&0\\0&-1\end{bmatrix}}$, The first column represents $\nabla _{x}c_{1}$ and second column represents $\nabla _{x}c_{2}$ .

Now we will calculate the values of $\nabla _{xx}c_{1}$ and $\nabla _{xx}c_{2}$ at the initial guess values.

$\nabla _{xx}c_{1}={\begin{bmatrix}0&1&0&0\\1&0&0&0\\0&0&0&0\\0&0&0&0\end{bmatrix}}$ , $\nabla _{xx}c_{2}={\begin{bmatrix}-2&0&0&0\\0&-2&0&0\\0&0&0&0\\0&0&0&0\end{bmatrix}}$

Now, $\nabla f(x)={\begin{bmatrix}x_{2}\\5+x_{1}\\0\\0\end{bmatrix}}={\begin{bmatrix}0\\5\\0\\0\end{bmatrix}}$. Similarly, $\nabla _{xx}f={\begin{bmatrix}0&1&0&0\\1&0&0&0\\0&0&0&0\\0&0&0&0\end{bmatrix}}$

In the above equation Matrix $'Z'$ and $'X'$ are basically diagonal matrices with the diagonal elements as follows:

$Z={\begin{bmatrix}z_{1}&0&0&0\\0&z_{2}&0&0\\0&0&z_{3}&0\\0&0&0&z_{4}\end{bmatrix}}$ and $X={\begin{bmatrix}x_{1}&0&0&0\\0&x_{2}&0&0\\0&0&x_{3}&0\\0&0&0&x_{4}\end{bmatrix}}$

As defined earlier, the matrix $'e'$ is simply a column vector of ones i.e. $e={\begin{bmatrix}1\\1\\1\\1\end{bmatrix}}$

Now since $W_{k}={\nabla _{xx}}^{2}L(x_{k},\lambda _{k},z_{k})={\nabla _{xx}}^{2}(f(x_{k})+c(x_{k})^{T}\lambda _{k}-z_{k})$ and we have chosen initial guess values of $\lambda _{1},\lambda _{2},z_{1},z_{2}as'0'$ , thus $W_{k}=\nabla _{xx}f$.

Now we have everything to substitute in the following equation and calculate the search directions and iterate till we fulfil the convergence criteria:

${\begin{bmatrix}W_{k}&\nabla c(x_{k})&-I\\\nabla c(x_{k})^{T}&0&0\\Z_{k}&0&X_{k}\end{bmatrix}}$${\begin{pmatrix}d_{k}^{x}\\d_{k}^{\lambda }\\d_{k}^{z}\end{pmatrix}}=$$-{\begin{pmatrix}\nabla f(x_{k})+\nabla c(x_{k})\lambda _{k}-z_{k}\\c(x_{k})\\X_{k}Z_{k}e-\mu _{j}e\end{pmatrix}}$

We can now implement BPOPT solver using MATLAB to get the following results as illustrated by :

Figure 2: Variation of x1 and x2 as iterations progress [9]

From figure 2 we see that the solver has converged to following optimal values of x1 and x2 in '13' iterations. $x1=4.3198$ and $x2=1.1575$

Figure 3: Variation of 'μ' as iterations progress [9]

In the figure 3 we can clearly see how the value of $'\mu '$ is decreasing as we progress towards the solution (increasing iterations).

Applications

IP methods have proven computationally to be viable alternative for the solution of several power engineering optimisation models. A variety of these are applied to a number of power systems models. 'Optimal power flow (OPF)' seeks the operation point of a power system that maximizes a given objective function and satisfies operational and physical constraints. Given real-time operating conditions and the large scale of the power system, it is demanding to develop algorithms that allow for OPF to be decomposed and efficiently solved on parallel computing systems . As compared to gradient based algorithms primal-dual based interior point methods have faster convergence rates.

A hydro-thermal coordination program deals with the problem of finding the scheduling and the production of every hydro and thermal plant of a system so that the customer demand is supplied at minimum cost with a certain level of security and all constraints related to the thermal and hydro subsystems are satisfied. A medium-term hydro-thermal coordination (MTHTC) problem results in a very large Non-Linear Mixed-Integer Programming problem IP optimisation codes are used to find the solution of these MTHTC problems.

IP methods are also used in Voltage collapse and reliability evaluation. Power flow unsolvability occurs when for a given set of active and reactive power loads the power flow equations have no real solution. This kind of problem may happen for instance when a heavily stressed system is subject to a severe contingency situation. The set of control actions in the IP algorithm includes rescheduling of active power of generators, adjustments on terminal voltage of generators, tap changes on LTC transformers, and as a last resort, minimum load shedding .

Fuel Planning and Multi-reservoir management are some other areas where interior point methods have proven to be useful. A long-term fuel planning problem aims to minimize the sum of fuel purchase and yearly storage costs, subject to meeting the generation requirements, the constraints associatedwith fuel supply contracts, and the constraints associated with storage of fuel.

Conclusion

Interior Point methods remain an active and fruitful area of research, although the frenetic pace that has characterized the area slowed down in recent years. Linear programming codes have become mainstream and continue to undergo development although they face stiff competition from the simplex method. Applications to quadratic programming show considerable promise, because of the superior ability of interior point approach to exploit problem structure efficiently. The influence on nonlinear programming theory and practice is yet to be determined even though substantial research has already been done on this topic . An extensive evaluation of line search options was conducted by and it was concluded that filter approach was the most robust approach. Many researchers are hopeful for the use of interior point approach in decomposition methods, though no rigorous comparative studies with alternative approaches have been performed. Many authors have also worked on a special subclass of non smooth optimization problems and a new interior point method has been proposed for this class. These methods extend a standard interior-point method for nonlinear and nonconvex problems by additional line-search and step length computation features and a problem-tailored stationarity test. A special type of generalized gradients is developed that is used at points of the nonsmooth problems that fail to be differentiable.