# Lecture 9: Methods using KKT conditions

## Sequential Quadratic Programming (SQP)

**Idea is to generate a sequence of quadratic optimization problems whose solutions approach the solution of the original problem**

Let us consider problem
$$
\min f(x)\\
\text{s.t. }h_k(x) = 0\text{ for all }k=1,\ldots,K,
$$
where the the objective function and the equality constraints are twice differentiable.

Because we know that the optimal solution of this problem satisfies the KKT conditions, we know that
$$
\left\{\begin{array}{l}
\nabla_xL(x,\lambda,\mu)=\nabla f(x) + \mu\nabla h(x) = 0\\
h(x) = 0
\end{array}\right.
$$
Let us assume that we have a current esimation for the solution of the equality constraints $(x^k,\mu^k)$, then according to the Newton's method for root finding (see e.g., https://en.wikipedia.org/wiki/Newton's_method), we have another solution $(x^k,\mu^k)^T+(p,v)^T$ of the problem by solving system of equations
$$
HL(x,\lambda,\mu)\left[\begin{align}p^T\\v^T\end{align}\right] = -\nabla L(x,\lambda,\mu).
$$
This can be written as
$$
\left[
\begin{array}{cc}
H_xL(x^k,\lambda,\mu^k)&\nabla h(x^k)\\
\nabla h(x^k)^T & 0
\end{array}
\right]
\left[\begin{array}{c}p^T\\v^T\end{array}\right] =
\left[
\begin{array}{c}
-\nabla_x L(x^k,\lambda,\mu^k)\\
-h(x^k)^T
\end{array}
\right].
$$

However, the above is just the solution of the quadratic problem with equality constraints
$$
\min \frac12 p^TH_xL(x^k,\lambda,\mu^k)p+\nabla_xL(x^k,\lambda,\mu^k)^Tp\\
\text{s.t. }h_j(x^k) + \nabla h_j(x^k)^Tp = 0. 
$$

## Intuitive interpretation

We are approximating the function quadratically around the current solution and the constraints are approximated linearly.

## Implementation

Define an optimization problem, where
* $f(x) = \|x\|^2$
* $h(x) = \sum_{i=1}^nx_i-n$


```python
def f_constrained(x):
    return sum([i**2 for i in x]),[],[sum(x)-len(x),x[0]**2+x[1]-2]
```


```python
print f_constrained([1,0,1])
print f_constrained([1,2,3,4])
```

    (2, [], [-1])
    (30, [], [6])



```python
import numpy as np
import ad



#if k=0, returns the gradient of lagrangian, if k=1, returns the hessian
def diff_L(f,x,m,k):
    #Define the lagrangian for given m and f
    L = lambda x_: f(x_)[0] + (np.matrix(f(x_)[2])*np.matrix(m).transpose())[0,0]
    return ad.gh(L)[k](x)
#Returns the gradients of the equality constraints
def grad_h(f,x):
    return  [ad.gh(lambda y:
                   f(y)[2][i])[0](x) for i in range(len(f(x)[2]))] 

#Solves the quadratic problem inside the SQP method
def solve_QP(f,x,m):
    left_side_first_row = np.concatenate((\
    np.matrix(diff_L(f,x,m,1)),\
    np.matrix(grad_h(f,x)).transpose()),axis=1)
    left_side_second_row = np.concatenate((\
    np.matrix(grad_h(f,x)),\
    np.matrix(np.zeros((len(f(x)[2]),len(f(x)[2]))))),axis=1)
    right_hand_side = np.concatenate((\
    -1*np.matrix(diff_L(f,x,m,0)).transpose(),
    -np.matrix(f(x)[2]).transpose()),axis = 0)
    left_hand_side = np.concatenate((\
                                    left_side_first_row,\
                                    left_side_second_row),axis = 0)
    temp = np.linalg.solve(left_hand_side,right_hand_side)
    return temp[:len(x)],temp[len(x):]exit
    
    

def SQP(f,start,precision):
    x = start
    m = np.ones(len(f(x)[2]))
    f_old = float('inf')
    f_new = f(x)[0]
    while abs(f_old-f_new)>precision:
        print x
        f_old = f_new
        (p,v) = solve_QP(f,x,m)
        x = x+np.array(p.transpose())[0]
        m = m+v
        f_new = f(x)[0]
    return x
```


```python
SQP(f_constrained,[0,0,0],0.0001)
```

    [0, 0, 0]
    [ 0.33333333  2.          0.66666667]
    [ 5.66666667 -1.66666667 -1.        ]
    [ 3.0199005  -0.11442786  0.09452736]
    [ 1.71499556  0.76156713  0.52343731]
    [ 1.12751769  1.0738341   0.79864821]
    [ 0.97581938  1.07078891  0.95339171]
    [ 1.01255467  0.97608252  1.01136281]
    [ 0.9939664   1.01237631  0.99365728]
    [ 1.00307036  0.99393273  1.00299691]





    array([ 0.99847908,  1.00306061,  0.99846031])



## Lagrangian methods -- "The original method of multipliers"

Let us again consider problem
$$
\min f(x)\\
\text{s.t. }h_k(x) = 0\text{ for all }k=1,\ldots,K,
$$
where the the objective function and the equality constraints are twice differentiable.

Define augmented Lagrangian function
$$
L_c(x,\mu) = f(x)+\mu h(x)+\frac12c\|h(x)\|^2.
$$
Above $c\in \mathbb R$ is a penalty parameter and $\mu \in \mathbb R^n$ is a multiplier.

Let us consider sequence of optimization problems
$$
\min_{x\in\mathbb R^n} f(x)+\mu_k h(x)+\frac{1}{2}c_k\|h(x)\|^2,
$$
where $c_{k+1}>c_k$ for $k=1,2,\ldots$.

Now, if $\mu_k=0$ for all $k=1,2,\ldots$, then we have a penalty function method, which solves the problem when $c_k\to \infty$.

However, it can be shown, that if we set $\mu_0$ randomly and keep on updating
$\mu_{k+1} = \mu_k-c_kh(x_k)$, then we can show that there exists $C>0$ such that of $c_k>C$, then the optimal solution of the augmented Langrangian solves the original problem!

### Example

Let us have optimization problem
$$
\min x_1^2+x_2^2\\
\text{s.t. }x_1+x_2-1=0.
$$

Now, the minimizatio of the augmented Lagrangian becomes
$$
\min_{x\in\mathbb R^n} x_1^2+x_2^2+\mu_k(x_1+x_2-1)+\frac12c_k(x_1+x_2-1)^2.\\
$$


```python
def f_constrained2(x):
    return sum([i**2 for i in x]),[],[sum(x)-1]
```


```python
def augmented_langrangian(f,x,mu,c):
    second_term = float(numpy.matrix(mu)*numpy.matrix(f(x)[2]).transpose())
    third_term = 0.5*c*numpy.linalg.norm(f(x)[2])**2
    return f(x)[0]-second_term+third_term
```


```python
from scipy.optimize import minimize
import numpy
def augmented_langrangian_method(f,start,mu0,c0):
    x_old = [float('inf')]*2
    x_new = start
    mu = mu0
    c = c0
    while numpy.linalg.norm(f(x_new)[2])>0.00001:
        res = minimize(lambda x:augmented_langrangian(f,x,mu,c),x_new)
        x_old = x_new
        mu = float(mu-numpy.matrix(c)*numpy.matrix(f(res.x)[2]).transpose())
        x_new = res.x
        c = 2*c
    return x_new,c
```


```python
from scipy.optimize import minimize
import numpy
def penalty_function_method(f,start,c0):
    x_old = [float('inf')]*2
    x_new = start
    c = c0
    while numpy.linalg.norm(f(x_new)[2])>0.00001:
        res = minimize(lambda x:augmented_langrangian(f,x,0,c),x_new)
        x_old = x_new
        x_new = res.x
        c = 2*c
    return x_new,c
```


```python
augmented_langrangian_method(f_constrained2,[0,0],1,1)
```




    (array([ 0.5,  0.5]), 2)




```python
penalty_function_method(f_constrained2,[0,0],1)
```




    (array([ 0.49999618,  0.49999618]), 262144)



## What is going on in here?

This is not completely trivial, unfortunately. If you want to read details, please see e.g., http://www.mit.edu/~dimitrib/Constrained-Opt.pdf.

599fb6285c21393c9a3aa8db0b1838dbecf2dee7
ipynb
Jupyter Notebook
Lecture 9, Methods using KKT conditions.ipynb
maeehart/TIES483
cce5c779aeb0ade5f959a2ed5cca982be5cf2316
CC-BY-3.0
2019-04-26T12:46:14.000Z
2021-11-23T03:38:59.000Z
Lecture 9, Methods using KKT conditions.ipynb
maeehart/TIES483
cce5c779aeb0ade5f959a2ed5cca982be5cf2316
CC-BY-3.0
Lecture 9, Methods using KKT conditions.ipynb
maeehart/TIES483
cce5c779aeb0ade5f959a2ed5cca982be5cf2316
CC-BY-3.0
2016-01-08T16:28:11.000Z
2021-04-10T05:18:10.000Z
Qwen/Qwen-72B
1. YES
2. YES