1
# -*- coding: utf-8 -*-
2
"""
3
Created on Tue Jun 21 14:11:29 2016
4

5
@author: Bossa
6

7
Algoritmo 6.3.2: Primal-Dual Path Following
8

9
Powered by Spider 2
10

11
"""
12
import numpy as np 
13

14
np.set_printoptions(precision=3)
15

16
class PrimalDualPF():
17
    '''
18
    Implementação do método Primal Dual Path Following para o problema 
19
    de programação linear 
20
        min  c*x
21
        s.a. Ax = b
22
             x >= 0   
23
    
24
    Parâmetros
25
    ----------
26
    A: matriz mxn com as variáveis de folga já inseridas.
27
    b: vetor nx1 
28
    c: vetor de custo mx1
29
    x: solução inicial do problema primal 
30
    w: solução inicial do problema dual
31
    s: folga dual associada a solução w
32
    '''
33
    def __init__(self, A, b, c, x, w, s):
34
        self.A = np.array(A) 
35
        # vetores serão arrays deitados
36
        self.b = np.array(b)
37
        # completando o vetor c com zeros
38
        self.c = np.zeros_like(A[0])
39
        self.c[0:len(c)] = np.array(c)
40
        self.x = np.array(x)
41
        self.w = np.array(w)
42
        self.s = np.array(s) 
43
        (m,n) = self.A.shape
44
        self.mu = np.dot(self.x,self.s)/n
45
        self.gamma = 1e-3
46
        
47
    def update(self):
48
        (m,n) = self.A.shape
49
        self.mu = np.dot(self.x,self.s)/n
50
        #Calculando a matriz Jacobiana
51
        J1 = np.concatenate( 
52
        (np.zeros((n,n)), np.transpose(self.A), np.eye(n)),axis=1)
53
        J2 = np.concatenate(
54
        (self.A, np.zeros((m,m)), np.zeros((m,n))),axis=1)
55
        J3 = np.concatenate(
56
        (np.diag(self.s), np.zeros((n,m)), np.diag(self.x)),axis=1)
57
        self.jacobian = np.concatenate((J1,J2,J3))
58
        # Calculando o lado direito, que será um array deitado
59
        D1 = self.c - self.A.T.dot(self.w) - self.s
60
        D2 = self.b - self.A.dot(self.x)
61
        D3 = -1*np.multiply(self.s,self.x)*np.ones(n) + \
62
               self.sigma*self.mu*np.ones(n)
63
        self.RHS = np.concatenate((D1,D2,D3)) 
64
        
65
    def step(self):
66
        (m,n) = self.A.shape 
67
        solucao = np.linalg.solve(self.jacobian,self.RHS) 
68
        # quebramsos a solucao do sistema linear em 
69
        # solucao = [px,pw,ps]
70
        px = solucao[0:n]
71
        pw = solucao[n:m+n]
72
        ps = solucao[m+n:m+2*n]
73
        alpha = 1.0
74
        backtracking = True
75
        while backtracking:
76
            # sejam 
77
            # (x_k+1,w_k+1,s_k+1) = (x_k,w_k,s_k) + alpha(px_k,pw_k,ps_k)
78
            novox = self.x + alpha*px
79
            novow = self.w + alpha*pw
80
            novos = self.s + alpha*ps 
81
            const = self.gamma*self.mu
82
            if validacao(novox,const) and validacao(novos,const):
83
                # se (x_k+1,s_k+1) > gamma*mu, tudo bem
84
                backtracking = False
85
            else:
86
                # senão, diminuímos o tamanho do passo
87
                alpha = 0.9*alpha
88
        self.x = novox
89
        self.w = novow
90
        self.s = novos
91
        
92
    def solve(self, TOL=1e-3,MAXSTEPS=1000, SIGMA=.2):
93
        '''Itera os passos de resolução até que a precisão TOL seja alcançada'''
94
        self.sigma = SIGMA
95
        self.caminho = [[],[]]
96
        for k in range(MAXSTEPS):
97
            if self.mu < TOL:
98
                break
99
            else:
100
                for j in range(2):
101
                    self.caminho[j].append(self.x[j])
102
                self.update()
103
                self.step()
104
            
105
    def solution(self):
106
        return (self.x,self.w,self.s)
107
            
108
def validacao(vector,constante):
109
    '''Verifica se cada entrada de um dado vetor é maior que a constante.'''
110
    test = True
111
    for x in vector:
112
        if x <= constante:
113
            test = False
114
            break
115
    return test
116
    
117
 
118

119