最优化方法（十五）有效集法解凸二次规划

#coding=gbk

'''
程序功能：
Date：2018.11.25
1、约束算法



'''

import sympy as sy
import numpy as np
import numpy.linalg as npl
import re
import matplotlib.pyplot as plt

#定义约束算法的类
class ConstrainedProgramming:
    

    #初始化
    def __init__(self,A,G,r,b,x):
        self.__A=A
        self.__G=G
        self.__r=r
        self.__b=b  #定义私有类变量
        

        self.__x=x
        n=len(self.__r)
        #print(n)
        
        self.__x=np.mat(self.__x).T  
        self.__indexM=self.__G.shape[0]+self.__A.T.shape[0]  #大矩阵的行
        self.__indexN=self.__G.shape[1]+self.__A.shape[1]  #大矩阵的列
        
        self.__calcResult()  

    #子函数计算处理
    def __calcResult(self):
    #fx=np.dot(fx,x)/2 # 
        fx=self.__objectiveFun()
        print('fx=',fx)
        

        
        AA=np.zeros((self.__indexM,self.__indexN))
        #print(AA)
        AA[0:self.__G.shape[0],0:self.__G.shape[1]]=self.__G
        AA[0:self.__G.shape[0],self.__G.shape[1]:self.__indexN]=self.__A 

        AA[self.__G.shape[0]:self.__indexM,0:self.__G.shape[1]]=self.__A.T
    
        #print(AA)
        bb=np.zeros((self.__indexM,1))
        #print(bb)
        bb[0:self.__G.shape[0]]=-self.__r

        bb[self.__G.shape[0]:self.__indexM]=self.__b
        #print(bb)
        self.__xx=npl.inv(AA).dot(bb) #
         
    
    
        
    #目标函数
    def __objectiveFun(self):
        fx=(self.__x.T).dot(self.__G).dot(self.__x)/2+self.__r.T.dot(self.__x)
        return fx
        #print(A)
        #print(b)
        
        
    #显示计算结果
    def show(self):
        print(self.__xx)
        print('---------------------------')
        print('无约束的解：\
',self.__xx[0:self.__G.shape[0]])
        print('无约束的乘因子：\
',self.__xx[self.__G.shape[0]:self.__indexM])

    #-------------------

#------------------------------
#定义约束算法有效集的类
class Youxiaoji:
    #初始化
    def __init__(self,x0,n,x):
        self.__n=n
        self.__x=x
        #print(self.__x)

       
        self.__num=10  #定义迭代次数
        self.__xinit=np.zeros((self.__num,2)) #计算的x点矩阵
        
        self.__xinit[0,:]=x0 #
        self.__d1fx=np.zeros((self.__num,2)) #一阶导数矩阵
        
        self.__calcResult()
    #------------------------------------------------------    
    #原函数的一阶导数   
    def __d1f(self): 
        f=self.__fun()
        y=[]
        for i in range(self.__n):
            y.append(sy.diff(f,self.__x[i]))
        y=np.array(y)
        print('d1f=',y)
        return y
    #------------------------------------------------------    
    #原函数
    def __fun(self):
        f=self.__x[0]**2+self.__x[1]**2-2*self.__x[0]-4*self.__x[1]
        return f
    #------------------------------------------------------
   #计算
    def __calcResult(self):
    
        f=self.__fun()
        print('f=',f)  
        d=self.__d1f()  
        for k in range(2): 
            self.__d1fx[0][k]=d[k].subs(self.__x[k],self.__xinit[0][k])
        print('----------------------')
        print('d1fx=\
',self.__d1fx)  #打印一阶导数值 



     #------------------------------------------------------   
    #显示结果
    def show(self):
       
        pass
     #------------------------------------------------------
 


#主函数
def main():
    
    #print(xx)
    #n=3
    #A=np.array([[1,2,-1],[1,-1,1]]).T
    #G=np.zeros((n,n))
    #r=np.zeros((n,1))
    #b=np.array([[4],[-2]])
    #x=sy.var('x1,x2,x3') #
    #for i in range(n):
    #    G[i][i]=2
    #--------------------------
    #n=2
    #A=np.array([[1,1]]).T
    #G=np.zeros((n,n))
    #r=np.zeros((n,1))
    #b=np.array([[1]])
    #x=sy.var('x1,x2') # 
    #for i in range(n):
    #    G[i][i]=2
    #c=ConstrainedProgramming(A,G,r,b,x) #定义类的对象
    #c.show()
    #---------约束初始条件-----------------
    n=2
    A=np.array([-1,0,0,-1]).reshape(2,2).T
    #print('A=',A)
    G=np.eye(n)*2
    #print('G=',G)
    r=np.array([-2,-4]).reshape(2,1)
    #print('r=',r)
    b=np.array([0,0]).reshape(2,1)
    #print('b=',b)
    x=sy.var('x1,x2') # 
    
    c=ConstrainedProgramming(A,G,r,b,x) #定义类的对象
    c.show()

    #--------------------------
    #有效集的程序初始化
    #x=sy.var('x1,x2') # 定义符号变量
    #n=len(x)
    #x0=np.array([[0,0]]) # 初始点
    #print('xinit=',x0)
    #v=Youxiaoji(x0,n,x)
    

main()