Решение тестовых задач с помощью Алгоритма 1

Краткое описание алгоритма

  1. Считывается матрица $A$ преобразования $\pi$, весовой вектор $w$. На основе последнего строится блочный порядок optimizationOrder.
  2. Строится кольцо многочленов $k[t,x]$ и идеал $I_t$ с помощью последовательной переработки столбцов матрицы $A$ в базисные многочлены.
  3. В базис $I_t$ добавляется нормирующий многочлен $t_0...t_d-1$.
  4. Ищется $gB$ -- базис Грёбнера идеала $I_t$.
  5. Из него выбираются элементы $gB_1=gB \cap k[x]$
  6. Достижимое решение $test$ записывается в виде монома, который затем сокращается базисом $gB_1$.

In [6]:
def TransformVectorsIntoPolynomials(vecs,orderr='lex'): #function for transforming vectors to Algorithm1 polynomials
    st=''
    variables=[]
    
    #Adding t>x
    for i in range(0,len(vecs[0])+1):
        st=st+'t'+str(i)+','
        variables.append('t'+str(i))
    
    #Adding x<t
    for i in range(0,len(vecs)-1):
        st=st+'x'+str(i+1)+','
        variables.append('x'+str(i+1))
    st=st+'x'+str(len(vecs))
    variables.append('x'+str(len(vecs)))
    
    #Polynomial Ring Construction 
    PR=PolynomialRing(QQ,st,order=orderr)
    pols=[]
    
    #Transformation from columns to Alg1 polynomials
    for i in range(0,len(vecs)):
        pol=''#result polynomial
        inds=set(range(0,len(vecs[i])))# unused indices
        first=False#is the variable the first in the monomial
        
        #generating minus part
        for j in range(0,len(vecs[i])):
            if(vecs[i][j]<0):
                if(not(first)):
                    first=True
                else:
                    pol=pol+'*'
                inds.remove(j)
                pol=pol+variables[j+1]#+1 because there is still t0
                
                if(vecs[i][j]<-1):
                    pol=pol+'^'+str(-vecs[i][j])
                    
        if(first):       #at least one negative is read
            pol=pol+'*'
        pol=pol+variables[len(vecs[0])+i+1]+'-'
        
        totalZero=True#if no positive entries
        first=False
        
        #generating plus part from unused indices
        for j in inds:
            if(vecs[i][j]!=0):
                if(not(first)):
                    first=True
                else:
                    pol=pol+'*'
                totalZero=False
                pol=pol+variables[j+1]#+1 because there is still t0
                if(vecs[i][j]>1):
                    pol=pol+'^'+str(vecs[i][j])
                    
        if(totalZero):
            pol=pol+'1'#if no positive entries
        
        pols.append(pol)#add to result set
        
        
    return [pols,PR]



def SolveOptimizationProblem(A,optimizationOrder,feasSol):
    
    [pols,PR]=TransformVectorsIntoPolynomials([A.column(i) for i in range(0,A.ncols())],optimizationOrder)
    #adding t0 polynomial
    pol=PR.gens()[0]
    for i in range(1,A.nrows()+1):
        pol=pol*PR.gens()[i]
    pol=pol-1

    pols.append(pol)
    print('******************The polynomial ring:')
    print(PR)    
    print('******************Generators of I_t:')
    for pol in pols:    
        print(pol)
    print('************************')
    
    #TermOrder('wdeglex',(1,1,1,1,
    #                     4,3,2,1,
    #                     7,5,3,1,
    #                    10,7,4,1)
    #pols

    I=PR.ideal(pols)

    #I
    gB=I.groebner_basis()

    print("Size of the Grobner basis of I_t: "+str(len(gB)))

    print("******************Valuable Groebner basis elements:")
    gB1=[]
    xVars=PR.gens()[(A.nrows()+1):]
    for g in gB:
        if(set(g.variables()) & set(xVars) == set(g.variables())):
            gB1.append(g)
            print(g)
    print("******************")

    
    polyn=1
    print(len(xVars),len(feasSol))
    for i in range(0,len(xVars)):
        polyn=polyn*xVars[i]^feasSol[i]

    print("Feasible Solution:"+str(feasSol))
    polynn=polyn.reduce(gB1)
    print("Optimal Solution:")
    print(polynn)
    
    return(str(polynn))

Транспортная задача [Sturmfels: ConvexPolytopes, стр. 41]

Решено верно.

Матрица $A$ и тестовые данные приведены ниже.

Задача:

$$ w^Tu \rightarrow min\\ s.t. Au=b, $$

где $w^T=[1,1,1,1,4,3,2,1,7,5,3,1,10,7,4,1],~b=[220,215,93,64,108,286,71,127]$.

Оптимальный план: $$ \begin{bmatrix} 108 & 112 & 0 & 0\\ 0 & 174 & 41 & 0\\ 0 & 0 & 30 & 63\\ 0 & 0 & 0 & 64 \end{bmatrix} $$


In [7]:
A=matrix([[1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0],
          [0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0],
          [0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0],
          [0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1],
          [1,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0],
          [0,1,0,0,0,1,0,0,0,1,0,0,0,1,0,0],
          [0,0,1,0,0,0,1,0,0,0,1,0,0,0,1,0],
          [0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1]])

    #first row is for t
optimizationOrder=TermOrder('lex',9)+TermOrder('wdeglex',(
                     1,1,1,1,
                     4,3,2,1,
                     7,5,3,1,
                    10,7,4,1))
feasibleSolution=vector([68,119,26,7,
        20,84,17,94,
        15,54,14,10,
        5,29,14,16])

SolveOptimizationProblem(A,optimizationOrder,feasibleSolution)


******************The polynomial ring:
Multivariate Polynomial Ring in t0, t1, t2, t3, t4, t5, t6, t7, t8, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x13, x14, x15, x16 over Rational Field
******************Generators of I_t:
x1-t1*t5
x2-t1*t6
x3-t1*t7
x4-t1*t8
x5-t2*t5
x6-t2*t6
x7-t2*t7
x8-t2*t8
x9-t3*t5
x10-t3*t6
x11-t3*t7
x12-t3*t8
x13-t4*t5
x14-t4*t6
x15-t4*t7
x16-t4*t8
t0*t1*t2*t3*t4*t5*t6*t7*t8 - 1
************************
Size of the Grobner basis of I_t: 127
******************Valuable Groebner basis elements:
x10*x13 - x9*x14
x6*x13 - x5*x14
x11*x13 - x9*x15
x7*x13 - x5*x15
x2*x13 - x1*x14
x3*x13 - x1*x15
x4*x13 - x1*x16
x8*x13 - x5*x16
x12*x13 - x9*x16
x6*x9 - x5*x10
x11*x14 - x10*x15
x7*x9 - x5*x11
x7*x14 - x6*x15
x2*x9 - x1*x10
x3*x9 - x1*x11
x3*x14 - x2*x15
x4*x9 - x1*x12
x4*x14 - x2*x16
x8*x9 - x5*x12
x8*x14 - x6*x16
x12*x14 - x10*x16
x7*x10 - x6*x11
x3*x10 - x2*x11
x4*x10 - x2*x12
x8*x10 - x6*x12
x2*x5 - x1*x6
x3*x5 - x1*x7
x4*x5 - x1*x8
x4*x15 - x3*x16
x8*x15 - x7*x16
x12*x15 - x11*x16
x3*x6 - x2*x7
x4*x6 - x2*x8
x4*x11 - x3*x12
x8*x11 - x7*x12
x4*x7 - x3*x8
******************
(16, 16)
Feasible Solution:(68, 119, 26, 7, 20, 84, 17, 94, 15, 54, 14, 10, 5, 29, 14, 16)
Optimal Solution:
x1^108*x2^112*x6^174*x7^41*x11^30*x12^63*x16^64
Out[7]:
'x1^108*x2^112*x6^174*x7^41*x11^30*x12^63*x16^64'

Просто тестовая задача с тетраэдром

Решено верно.

Минимизировать $w^Tx$ при $w \in R^3,~ b \in Z^3, ~x \in N^3 $ и ограничениях $A x <= b, ~x>=0$.

$$A=\begin{bmatrix}1 & 1 & 1\\ -2 & -1 & 0\\ 0 & -1 & 0 \end{bmatrix}\\ ~b=[12,-6,-6]^T,\\ w=[1,1,1]^T$$

.

Преобразуем ограничения: $Ax +s = b, s>=0$.

Достижимое решение (с учётом $s$): $[2,10,0,0,16,4]^T$.


In [8]:
A=matrix([[1,1,1,1,0,0],[-2,-1,0,0,1,0],[0,-1,0,0,0,1]])
#initial weights 'wdeglex',(1,5,4)
#lagrange: weight of the constraint is -1
#+constraints to w??
dw=vector(2*A.row(0)+A.row(1)+A.row(2))
w=tuple(vector([1,1,1,0,0,0])+dw)
optimizationOrder=TermOrder('lex',4)+TermOrder('wdeglex',w)
feasibleSolution=vector([2,10,0,0,2,4])
SolveOptimizationProblem(A,optimizationOrder,feasibleSolution)


******************The polynomial ring:
Multivariate Polynomial Ring in t0, t1, t2, t3, x1, x2, x3, x4, x5, x6 over Rational Field
******************Generators of I_t:
t2^2*x1-t1
t2*t3*x2-t1
x3-t1
x4-t1
x5-t2
x6-t3
t0*t1*t2*t3 - 1
************************
Size of the Grobner basis of I_t: 10
******************Valuable Groebner basis elements:
x2^2*x6^2 - x1*x4
x2*x5*x6 - x4
x3 - x4
x1*x5 - x2*x6
******************
(6, 6)
Feasible Solution:(2, 10, 0, 0, 2, 4)
Optimal Solution:
x1^3*x2^6*x4^3
Out[8]:
'x1^3*x2^6*x4^3'

Bounded Knapsack Problem (BKP)

Решение верное(проверено графически в MATLAB).

Формулировка взята отсюда: http://www.or.deis.unibo.it/kp/Chapter3.pdf (целая книга про задачи о рюкзаке).

Максимизировать $w^Tx$ при $w \in R^3,~ b \in Z^4, ~x \in N^3 $ и ограничениях $A x <= b, ~x>=0$.

$$A=\begin{bmatrix} 1 & 3 & 5\\ 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end{bmatrix}\\ ~b=[10,6,4,2]^T,\\ w=[10,15,11]^T. $$

Преобразуем ограничения: $Ax +s = b, s>=0$.

Преобразованный $-w$ (задача $min -w^Tx$): $$-w+10Constraint_1+Constraint_2+Constraint_3+Constraint_4=[2,19,45,11,1,1,1].$$

Достижимое решение (с учётом $s$): $[0,3,0,1,6,1,2]^T$.


In [9]:
A=matrix([[1,3,5,1,0,0,0],[1,0,0,0,1,0,0],[0,1,0,0,0,1,0],[0,0,1,0,0,0,1]])
#initial weights 'wdeglex',(1,5,4)
#lagrange: weight of the constraint is -1
#+constraints to w??
dw=vector(10*A.row(0)+A.row(1)+A.row(2)+A.row(3))
w=tuple(-vector([10,15,11,0,0,0,0])+dw)
optimizationOrder=TermOrder('lex',5)+TermOrder('wdeglex',w)
feasibleSolution=vector([0,3,0,1,6,1,2])
SolveOptimizationProblem(A,optimizationOrder,feasibleSolution)


******************The polynomial ring:
Multivariate Polynomial Ring in t0, t1, t2, t3, t4, x1, x2, x3, x4, x5, x6, x7 over Rational Field
******************Generators of I_t:
x1-t1*t2
x2-t1^3*t3
x3-t1^5*t4
x4-t1
x5-t2
x6-t3
x7-t4
t0*t1*t2*t3*t4 - 1
************************
Size of the Grobner basis of I_t: 17
******************Valuable Groebner basis elements:
x4^5*x7 - x3
x1*x4^4*x7 - x3*x5
x3*x5^2 - x1^2*x4^3*x7
x3*x6 - x2*x4^2*x7
x4^3*x6 - x2
x1*x4^2*x6 - x2*x5
x2*x5^2 - x1^2*x4*x6
x4*x5 - x1
******************
(7, 7)
Feasible Solution:(0, 3, 0, 1, 6, 1, 2)
Optimal Solution:
x1^6*x2*x4*x6^3*x7^2
Out[9]:
'x1^6*x2*x4*x6^3*x7^2'

In [ ]: