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))
Решено верно.
Матрица $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)
Out[7]:
Решено верно.
Минимизировать $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)
Out[8]:
Решение верное(проверено графически в 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)
Out[9]:
In [ ]: