数理計画における並列計算 (経営科学のニューフロンティア) 　朝倉書店

```{r}

#数理計画における数理最適

#2次輸送問題 p.47

library(dplyr)

m=10;

n=m;

vec=sample(1:100,m,replace=T)

X=vec

sita=sample(10:20,m)

C=sample(10:20,m)

vec=sample(1:100,m,replace=T)

A=t(t(vec))%*%vec

b=apply(A,1,sum)+100000;z=rep(10,m)


w=0.01

ite=1000

for(l in 1:ite){
  
X_sub=X;  
  
for(i in 1:length(X)){
  
X_sub[i]=-(sum(A[i,]*z)+C[i])/sita[i]  
  
} 


X=X_sub;    
  
z_pre=z;  
  
for(i in 1:length(z)){
  
z[i]=max(c(0,z[i]-w*(b[i]-sum(A[i,]*X))/sum(A[i,]^2/sita)))  
  
}  


  
cost=sum(sita*X^2)/2+sum(C*X)

print(cost)
  
}

A%*%X<=b


```

```{r}

#2次輸送問題 p.51

library(dplyr)

m=10;

n=m;

vec=sample(1:100,m,replace=T)

X=t(t(vec))%*%vec+1000

sita=array(0,dim=c(m,n))

C=array(0,dim=c(m,n))

for(i in 1:m){
for(j in 1:n){  
  
sita[i,j]=sample(10:20,1)

C[i,j]=sample(20:30,1)
    
}
}


w=0.001

u=rep(1,m);v=rep(2,n)

w_mat=array(2,dim=c(m,n))

ite=1000

for(l in 1:ite){
  
w_sub=w_mat;u_sub=u;v_sub=v;X_sub=X  
  
alpha=apply(X,1,sum);beta=apply(X,2,sum)
  
for(i in 1:m){
for(j in 1:n){  
 
X[i,j]=(w_mat[i,j]-u[i]-v[j]-C[i,j])/sita[i,j]   
  
}
}

u_pre=u
  
for(i in 1:m){
  
u[i]=u[i]-w*(alpha[i]-apply(X,1,sum)[i])/(apply(1/sita,1,sum)[i])  
  
}  

X2=X;

for(i in 1:m){
  
X2[i,]=X[i,]-(u[i]-u_pre[i])/sita[i,]  
  
}

v_pre=v

for(j in 1:n){
  
v[j]=v[j]-w*(beta[j]-apply(X2,2,sum)[j])/(apply(1/sita,2,sum)[j])  
  
}

X3=X2;

for(j in 1:n){
  
X3[,j]=X2[,j]-(v[j]-v_pre[j])/sita[,j]  
  
}

for(i in 1:m){
for(j in 1:n){  
  
w_mat[i,j]=max(c(0,w_mat[i,j]-w*sita[i,j]*X3[i,j])) 
   
}
}    
  
  
#print(sum(abs(w_mat-w_sub)+abs(u_sub-u)+abs(v_sub-v)+abs(X_sub-X)))

print(sum(sita*X^2)/2+sum(C*X))

}

sum(alpha)==sum(beta)

```
```{r}

#PVT method p.114

alpha=0.01

m=5

p=sample(10:20,m,replace=T)


f=function(x){
  
return(sum(abs(x)^3))  
  
}

ite=10000

X=sample(1:10,m,replace=T)

h=0.001;B=diag(1,m)

g_pre=rep(0,m);g=rep(0,m)

cost=p%*%B%*%t(t(p))/2+g%*%t(t(p))+f(X)

cost=as.numeric(cost);cost_pre=cost+10

for(l in 1:ite){
 
if(abs(cost_pre-cost)>0.01){  
  
cost_pre=cost  
   
g_pre=g  

X_pre=X
  
g=c()

vec=X

for(i in 1:length(X)){
 
vec_sub=vec;vec_sub[i]=vec_sub[i]+h  
   
g=c(g,(f(vec_sub)-f(vec))/h)  
  
}

H=array(0,dim=c(m,m))

for(i in 1:m){
for(j in 1:m){  

vec_sub1=vec;vec_sub2=vec;vec_sub1[i]=vec_sub1[i]+h;vec_sub2[j]=vec_sub2[j]+h 

vec_sub3=vec;vec_sub3[i]=vec_sub3[i]+h;vec_sub3[j]=vec_sub3[j]+h
   
H[i,j]=(f(vec_sub3)-f(vec_sub1)-f(vec_sub2)-f(vec))/(h^2)   
  
}
}  

X=X+alpha*(-solve(H)%*%g)  

S=c(X-X_pre)

y=c(g-g_pre)

mat=-(S%*%B%*%t(t(S)))/alpha

#u=y/as.numeric(t(y)%*%t(t(S)))+g*sqrt(-alpha/as.numeric(mat%*%(t(y)%*%t(t(S)))))

B=B+t(t(y))%*%t(y)/as.numeric((t(y))%*%t(t(S)))-B%*%t(t(S))%*%t((S))%*%B/as.numeric(S%*%B%*%t(t(S)))
 
cost=p%*%B%*%t(t(p))/2+g%*%t(t(p))+f(X)

cost=as.numeric(cost)

print(cost)

}
  
}  


```

```{r}

#p.165

library(dplyr)

m=5;eta=0.001

C1=sample(1:10,m,replace=T);C2=sample(1:10,m,replace=T);

b1=sample(5:10,m,replace = T)

b2=sample(5:10,m,replace = T)

b0=sample(500:1000,m,replace = T)

A=matrix(sample(0:1,m*m,replace=T),ncol=m)

ite=10000

X=rep(1,m)

vec1=sample(10:100,m,replace=T);vec2=sample(10:100,m,replace=T)

sita1=t(t(vec1))%*%vec1/1000;sita2=t(t(vec2))%*%vec2/1000

s=rep(1,m);w1=vec1;w2=vec2;x1=vec2;x2=vec2;y=rep(1,m);z1=vec2;z2=vec2

mu_val=sort(seq(0.1,10,0.1))

cost_vec=c()

cost_pre=10^10+10;cost=10^10

v_pre=10+1;v=10

mu=array(100,dim=c(m,1));v=10

for(l in 1:ite){
  
if(abs(v-v_pre)>0.01){  
  
v_pre=v  
  
s_pre=s;w1_pre=w1;w2_pre=w2;x1_pre=x1;x2_pre=x2;y_pre=y;z1_pre=z1;z2_pre=z2  
  
D1=solve(sita1%*%diag(x1)+diag(w1))%*%diag(x1)

D2=solve(sita2%*%diag(x2)+diag(w2))%*%diag(x2)  

e1=solve(sita1%*%diag(x1)+diag(w1))%*%(mu-diag(x1)%*%w1)

e2=solve(sita2%*%diag(x2)+diag(w2))%*%(mu-diag(x2)%*%w2)

AD1A=solve(A%*%D1%*%t(A)+diag(10^(-10),m))

AD2A=solve(A%*%D2%*%t(A)+diag(10^(-10),m))

G=diag(y)%*%((D1-D1%*%t(A)%*%(AD1A)%*%A%*%D1)+(D2-D2%*%t(A)%*%(AD2A)%*%A%*%D2))+diag(s)

hs=diag(y)%*%((e1-D1%*%t(A)%*%(AD1A)%*%A%*%e1)+(e2-D2%*%t(A)%*%(AD2A)%*%A%*%e2))+mu-diag(s)%*%y

dy=solve(G,hs);ds=solve(diag(y))%*%(mu-diag(s)%*%c(y+dy))  

dz1=(AD1A)%*%A%*%(e1-D1%*%dy)

dz2=(AD2A)%*%A%*%(e2-D2%*%dy)

dx1=e1-D1%*%(t(A)%*%dz1+dy);dx2=e2-D2%*%(t(A)%*%dz2+dy)

dw1=solve(diag(x1))%*%(mu-diag(w1)%*%c(x1+dx1))

dw2=solve(diag(x2))%*%(mu-diag(w2)%*%c(x2+dx2))

s=c(s+eta*ds)

y=c(y+eta*dy)

mu=diag(s)%*%y

x1=c(x1+eta*dx1);x2=c(x2+eta*dx2)

w1=c(w1+eta*dw1);w2=c(w2+eta*dw2)

z1=c(z1+eta*dz1);z2=c(z2+eta*dz2)

v0=sum(abs(mu-diag(s)%*%y));v1=sum(abs(mu-diag(x1)%*%w1));v2=sum(abs(mu-diag(x2)%*%w2))

v=max(c(v0,v1,v2))

#cost=sum(abs(x1-x1_pre))+sum(abs(x2-x2_pre))+sum(abs(w1_pre-w1))+sum(abs(w2_pre-w2))+sum(abs(z1_pre-z1))+sum(abs(z2_pre-z2))+sum(abs(s_pre-s))+sum(abs(y_pre-y))

cost=as.numeric(t(x1)%*%sita1%*%t(t(x1)))/2+sum((C1)*x1)+as.numeric(t(x2)%*%sita2%*%t(t(x2)))/2+sum((C2)*x2)

print(paste0("cost:",cost," v:",v))

cost_vec=c(cost_vec,cost)

}

}
  
plot(cost_vec)


```
```{r}

#数理計画入門

#BFKG method p.113

library(dplyr)

f=function(x){
  
  return(sum(c(exp(abs(x[1]-2))-1,exp(abs(x[2]-1))-1)))
  
}


B=diag(1,2)

x=rep(3,2)

ite=2000

alpha=0.001

beta=-1

h=0.01

for(l in 1:ite){
  
df=c();vec=x

for(i in 1:2){
  
vec_sub=vec;vec_sub[i]=vec_sub[i]+h  
  
df=c(df,sum((f(vec_sub)-f(vec))/h))  
  
}
  
d=-solve(B)%*%df

x=x+alpha*d

y=c()

x_pre=vec

vec_pre=x_pre;vec=x

for(i in 1:2){
  
vec_pre_sub=vec_pre;vec_sub=vec

vec_pre_sub[i]=vec_pre_sub[i]+h;vec_sub[i]=vec_sub[i]+h
  
y=c(y,(f(vec_pre_sub)-f(vec_pre))/h-(f(vec_sub)-f(vec))/h)  
  
}

s=c(x-vec)

beta=as.numeric(t(y)%*%t(t(s)))

if(beta>0){

B=B+t(t(y))%*%t(y)/as.numeric(t((y))%*%t(t(s)))-B%*%t(t(s))%*%t(s)%*%B/as.numeric(t(s)%*%B%*%t(t(s)))
  
}

print(x)

}

```
```{r}

#p.149 lambdaについて疑問

m=100;n=60

A=matrix(sample(1:100,m*n,replace=T),ncol=m)/10000

b=sample(1:10,n,replace=T)*100

sita=sample(1:100,m,replace=T)/1000;C=-sample(1:20,m,replace=T)*10

lambda=10^(-6)



ite=20000

y=array(1,dim=dim(A))

z=array(1,dim=dim(A))

x=solve(t(A)%*%A+diag(1,m))%*%t(A)%*%b;cost=-10^20

r=rep(1,n);u=rep(1,n)

guzai=apply(A^2,1,sum)


for(l in 1:ite){

if(cost<0){  

cost_pre=cost
  
for(j in 1:length(x)){

x[j]=max(c(0,(sum(lambda*y[,j]-z[,j])-C[j])/(2*lambda+sita[j])))  

#x[j]=max(c(0,(2*lambda*x[j]-C[j]-sum((u+lambda*r/guzai)*A[,j]))/(2*lambda+sita[j])))  

}

r=A%*%x-b


u=(lambda*r+apply(A*z,1,sum))/guzai

#y=t(t((z-A*c(u))/lambda)+c(x))

y=t(t(-A*c(r)/c(guzai))+c(x))

z=z-t(t(y)-c(x))*lambda

cost=as.numeric(sum(sita*x^2)/2+sum(C*x)+sum(-t(t(y)-c(x))*z)+lambda*sum(t(t(y)-c(x))^2)/2)

print(cost)    

}

}  

```

```{r}

#p.198

m=5

vec=sample(1:10,m,replace=T)

sita=t(t(vec))%*%t(vec)

C=sample(1:20,m,replace = T)

val=sample(1:10,m,replace=T)/10

f=function(X){
  
X=c(X)  
  
return(sum(diag(sita)*X^2)+sum((sita-diag(diag(sita)))*(t(t(X))%*%t(X)))+sum(C*(X-val))/2)  
  
}


d_mat=array(10,dim=c(m,m-1))

x_mat=array(10,dim=c(m,m))

ite1=10;ite2=300000;ite3=300000

eta=0.00001

w=rep(1,m+1)

h=0.01

x=vec

for(l in 1:ite1){
  
for(i in 1:m){
  
no=c(1:m);no=no[-i]

vec=x;x_sub=x

cost=1000

for(s in 1:ite2){
  
if(abs(cost)>1){  
  
vector=x_sub  
  
for(j in no){
  
vec_sub1=vector;vec_sub1[j]=vec_sub1[j]+h  

vec_sub2=vector;vec_sub2[j]=vec_sub2[j]-h
  
x_sub[j]=x_sub[j]-((f(vec_sub1)-f(vector))/h)*eta
  
}
  
cost=f(x_sub)

if(cost<0){x_sub=vector}
  
}  
  
}  

print(paste0("no.",i,",",as.numeric(f(x_sub))))
  
d=x_sub-x;d_mat[i,]=d[no];x_mat[i,]=x   

}  
  
cost=1000; 
  
eta2=10^(-8)

for(s in 1:ite3){
  
if(abs(cost)>1){  
  
vec=w 
 
for(j in 1:(m+1)){

vec_sub=vec;vec_sub[j]=vec_sub[j]+h  

vec_sub_sub=vec;vec_sub_sub[j]=vec_sub_sub[j]-h

x1=apply((rbind(x,x_mat))*c(vec_sub),2,sum)

x2=apply((rbind(x,x_mat))*c(vec),2,sum)

x3=apply((rbind(x,x_mat))*c(vec_sub_sub),2,sum)

df=(f(x1)-f(x2))/h 

ddf=(f(x1)-2*f(x2)+f(x3))/(h^2)

w[j]=w[j]-df*eta2 
  
}

cost=f(apply((rbind(x,x_mat))*c(w),2,sum))

if(cost<0){w=vec}

cost=f(apply((rbind(x,x_mat))*c(w),2,sum))

print(cost)

}

} 

print(f(apply((rbind(x,x_mat))*c(w),2,sum)))

x=apply((rbind(x,x_mat))*c(w),2,sum)  
  
}


```