線形計画法を用いて当直の最適な割り当てを考える

R 統計

MikuHatsune2013-02-07

シンプレックス法やIsotonic regressionなど、線形計画法の応用をいろいろやっているが、Kashiwa.R#5でやったネタをみた先生から、

最近、内科全体の当直表が○科で重なりまくり。
文句を言ったら、重ならないようには作れないだと!どうみても余裕がある。

との相談を受けたので、プログラミングでできるようにやってみる。
　
n日の当直d_j;j=1,2,\dots,nをm人の医師p_i;i=1,2,\dots,mでカバーする。
ある当直日d_jには、医師がL_{d_j}人待機していないといけない。
ある医師p_iは、n日の間で最低w_{p_i}回、最大でW_{p_i}回当直をしなければならない。
ある医師p_iは所属診療科D_{p_i}の属性がある。
ある医師p_iがd_jに当直をしたら、d_0の中日を設けて休みをとらなければならない。
　
という設定でやろう。使うのは、lpSolveパッケージの

lp.assign

でやりたかったけど、複数の当直を割り当てるのがうまくいかなかったので、線形計画法の解を整数指定することでやろう。

lp(int.vec = ...)

でできる。
　
lpに渡すときには、行列が長々とベクトル化されて計算されるようなので、それを踏まえて制限行列を作成する。
医師列p_iと当直日行d_jの各項に、割り当て項x_{i,j}=\{0,1\}をいれる。
\begin{matrix}&d_1&d_2&\dots&d_j&\dots&d_n\\p_1&x_{1,1}&x_{1,2}&\dots&x_{1,j}&\dots&x_{1,n}\\p_2&x_{2,1}&x_{2,2}&\dots&x_{2,j}&\dots&x_{2,n}\\\vdots&&&&&&\\p_i&x_{i,1}&x_{i,2}&\dots&x_{i,j}&\dots&x_{i,n}\\\vdots&&&&&&\\p_m&x_{m,1}&x_{m,2}&\dots&x_{m,j}&\dots&x_{m,n}\end{matrix}
これをbyrowでベクトルに崩して長い制限ベクトルを作ったとする。
制限行列は(+1を省略)
\begin{matrix}w_{p_1}\leq&x_{1,1}&x_{1,2}&\dots&x_{1,j}&\dots&x_{1,n}&&&&&&&&&&&&&&&&&&&&&&&&&&&&&\leq W_{p_1}\\w_{p_2}\leq&&&&&&&x_{2,1}&x_{2,2}&\dots&x_{2,j}&\dots&x_{2,n}&&&&&&&&&&&&&&&&&&&&&&&\leq W_{p_2}\\\vdots&&&&&&&&&&&&&&\ddots&&&&&&&&&&&&&&&&&&&&&\vdots\\w_{p_i}\leq&&&&&&&&&&&&&&&&&&x_{i,1}&x_{i,2}&\dots&x_{i,j}&\dots&x_{i,n}&&&&&&&&&&&&\leq W_{p_i}\\\vdots&&&&&&&&&&&&&&&&&&&&&&&&&\ddots&&&&&&&&&&\vdots\\w_{p_m}\leq&&&&&&&&&&&&&&&&&&&&&&&&&&&&&x_{m,1}&x_{m,2}&\dots&x_{m,j}&\dots&x_{m,n}&\leq W_{p_m}\\&x_{1,1}&&&&&&x_{2,1}&&&&&&\dots&&&&&x_{i,1}&&&&&&\dots&&&&&x_{m,1}&&&&&&\leq L_{d_1}\\&&x_{1,2}&&&&&&x_{2,2}&&&&&\dots&&&&&&x_{i,2}&&&&&\dots&&&&&&x_{m,2}&&&&&\leq L_{d_2}\\&\vdots&&\ddots&&&&&&\ddots&&&&&&&&&&&\ddots&&&&&&&&&&&\ddots&&&&\vdots\\&&&&x_{1,j}&&&&&&x_{2,j}&&&&&\dots&&&&&&x_{i,j}&&&&&\dots&&&&&&x_{m,j}&&&\leq L_{d_j}\\&\vdots&&&&\ddots&&&&&&\ddots&&&&&&&&&&&\ddots&&&&&&&&&&&\ddots&&\vdots\\&&&&&&x_{1,n}&&&&&&x_{2,n}&&&&&\dots&&&&&&x_{i,n}&&&&&\dots&&&&&&x_{m,n}&\leq L_{d_n}\\0\leq&x_{1,1}&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&\leq 1\\0\leq&&x_{1,2}&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&\leq 1\\\vdots&&&\ddots&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&\vdots\\0\leq&&&&&&&&&&&&&&&&&&&&&x_{i,j}&&&&&&&&&&&&&&\leq 1\\\vdots&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&\ddots&\vdots\\0\leq&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&x_{m,n}\leq 1\end{matrix}
　
という感じでやると、こんな感じに割り当てられる。
中日は2日にした。
このやり方で問題なのが、当直をやります、と希望した前後の日は、かなり大きな負のスコアを入れることで、割り当てられないようにしていることだ。
この一週間のうち、4日間は適当に割り当ててくれたらいいよ、というような要望を強制的に変換してしまうので、線形計画法で解をd_0日とるような解法を組み込みたい。
今回は希望調査行列から計算されて出てきた割り当て行列について、中日d_0、同一診療科医師の当直D_{p_a}=D_{p_b};a\not= bを逐一チェックして、ダメなときは再計算する方法にした。
これだと、最適な割り当てが一生見つからない（そもそも制限がない状況で最適解がないときもあるが）ことになりうるので、まあどうにかしよう。

time0 <- as.Date("2013-11-01")
time1 <- seq(time0, len=5, by="1 month")[2]
seq(time0, time1, by = 1)
days <- head(seq(time0, time1, by = 1), -1L)
holidays <- as.integer(weekdays(days) == "日曜日")

mat <- read.delim("clipboard", header=TRUE)
dept <- mat[, 1]
doctors <- nrow(mat)
mat <- as.matrix(mat[, -1])
library(lpSolve)
#パラメータ群
doctors <- 20 #医師の人数
n_department <- 12 #診療科の数
n_toutyoku <- 4 #医師の最大当直回数(月)
n_doctors <- 2 #一日あたりの必要当直医師数
min_toutyoku <- 1 #医師の最小当直回数(月)
nakabi0 <- 2 #一回当直したあと、次の当直までの中日
penalty <- -100 #中日に当直が当たらないようにするペナルティスコア

dept <- sort(sample(seq(n_department), size=doctors, replace=TRUE)) #医師の所属科
mat <- matrix(0, nr = doctors, nc = length(days)) #希望日程行列
for(i in seq(doctors)){
	mat[i, ] <- rpois(length(days), 1) * sample(c(-1, 0, 1), size=length(days), prob=c(0.05, 0.9, 0.05), replace=TRUE)
}
#当直前後の日にあらかじめペナルティをいれておく
for(i in seq(doctors)){
	tmpv <- c(rep(0, nakabi0), mat[i, ], rep(0, nakabi0))
	for(j in seq(length(days)) + nakabi0){
		if(tmpv[j] > 0){
			tmpv[c((j - nakabi0):(j - 1), (j + 1):(j  + nakabi0))] <- penalty
		}
	}
	mat[i, ] <- tmpv[seq(length(days)) + nakabi0]
}

#制限行列の作成
Q1 <- matrix(c(mapply(rep, diag(1, doctors), len=length(days))), nr=doctors, byrow=TRUE) #各医師の当直日数
Q2 <- t(apply(diag(1, length(days)), 1, rep, doctors)) #一日に必要な当直医師数
Q3 <- diag(1, prod(dim(mat))) #各項の割当は0か1

constQ <- rbind(Q1, Q2, Q3, Q1)
duty_doc <- rep(n_toutyoku, length=doctors) #医師の当直回数。Q1に相当
duty_days <- rep(n_doctors, length=length(days)) #一日に拘束される医師数。Q2に相当
min_vec <- rep(min_toutyoku, doctors) #各医師の当直日数の最低数
const_vec <- c(duty_doc, duty_days, rep(1, nrow(Q3)), min_vec)

#制限行列の不等式
direc <- unlist(mapply(rep, c("<=", "==", "<=", ">="), c(doctors, length(days), nrow(Q3), doctors)))
#cbind(const_vec, direc) #確認

#中日が何日以上か空いていないと採用しない関数
nakabi <- function(vec, days = 4){
	v0 <- rev(which(vec != 0))
	if(!all(v0 == 0)){
	vdiff <- head(v0, -1L) - tail(v0, -1L)
	return(all(vdiff > days))
	}else{
		return(FALSE)
	}
}

#同一診療科の医師が同日に当直に入っていたら採用しない関数
kaburi <- function(assign_vec, department){
	dvec <- department[which(assign_vec == 1)]
	return(length(dvec) == length(unique(dvec))) #被りはFALSE
}

#探しにかかる
for(i in seq(10000)){
	id <- sample(seq(nrow(mat)))
	mat0 <- mat[id, ]
	res <- lp("max", c(t(mat0)), constQ, direc, const_vec, int.vec=seq(prod(dim(mat))))
	duty0 <- matrix(res$solution, nr=doctors, byrow=TRUE)[id, ]
	#if(all(apply(duty0, 1, nakabi, days = nakabi0))){ #中日だけチェックする
	#if(all(apply(duty0, 2, kaburi, department = dept))){ #診療科の被りだけチェックする
	#ゴリゴリプロットする
	par(mar=c(3, 6, 2, 2), xpd=TRUE, mgp=c(0.1, 1.5, 0))
	image(duty0, xlab="", ylab="", axes=FALSE, col=terrain.colors(12))
	#所属診療科を色で表示する
	xl <- seq(par()$usr[1], par()$usr[2], length=doctors + 1)
	yl <- c(par()$usr[3], par()$usr[3] - 1/max(doctors))
	for(i in seq(dept)){
		polygon(xl[c(i:(i+1), (i+1):i)], yl[rep(1:2, each=2)], col=rainbow(max(dept))[dept[i]], border=FALSE, xpd=TRUE)
	}
	#日曜日に目印をつける
	xl <- c(par()$usr[1], par()$usr[1] - 1/max(doctors))
	yl <- seq(par()$usr[3], par()$usr[4], length=length(days) + 1)
	for(i in seq(length(days))){
		polygon(xl[rep(1:2, each=2)], yl[c(i:(i+1), (i+1):i)], col=holidays[i], border=FALSE, xpd=TRUE)
	}
	axis(1, at=seq(0, 1, length=doctors), seq(doctors), las=2, tck=0, fg=NA)
	axis(2, at=seq(0, 1, length=length(days)), as.character(days), las=1, tck=0, fg=NA)
	
	print(i)
	print(colSums(duty0)) #一日の必要当直人数の確認
	print(rowSums(duty0)) #各医師の当直回数の確認	
	if(all(apply(duty0, 1, nakabi, days = nakabi0)) & all(apply(duty0, 2, kaburi, department = dept))){
		break #いいのがあったらやめる
	}
}

いい感じになったときの希望調査

dept	V1	V2	V3	V4	V5	V6	V7	V8	V9	V10	V11	V12	V13	V14	V15	V16	V17	V18	V19	V20	V21	V22	V23	V24	V25	V26	V27	V28	V29	V30
1	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	-1	0	0	2
1	0	0	0	0	0	-1	0	0	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0
1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
2	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	-1	0	0	0	0	0	0	0	0	0	0	0
2	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
3	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	0	-1	0	0	0	0	0	0
5	0	0	0	0	0	0	0	0	0	0	2	0	0	0	0	0	0	1	0	0	0	0	1	0	0	0	0	0	0	0
5	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
7	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	-2	0	0	0	0	0	0	0	0	0	0	0	0	0
8	0	-2	0	0	3	0	0	0	0	0	0	0	0	0	2	0	0	0	0	0	1	0	0	0	0	0	0	-2	0	0
8	0	0	0	0	0	-1	0	0	0	0	0	0	0	0	0	0	0	0	3	0	0	0	0	0	0	0	0	0	0	0
8	0	0	0	0	0	0	0	0	0	0	0	0	0	2	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
9	0	0	0	0	0	2	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1
9	0	0	0	0	0	0	2	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	2	0	0	0	0
10	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	-1	0	0	0	0	0	0	0	0	0	0	0	0	0	0
10	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	1	0	0
11	0	0	2	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
11	0	0	0	0	0	0	0	0	0	0	0	0	2	0	0	0	0	0	0	0	-1	0	0	0	0	0	0	0	-1	0
12	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	-2	0	0	0	0	0	0	0	0	0	0	1	0	0	0
12	0	0	-1	0	0	-2	0	0	0	0	0	0	-2	0	0	0	0	0	-1	-1	0	0	0	0	0	0	3	0	0	0

そのときの当直割り当て

V1	V2	V3	V4	V5	V6	V7	V8	V9	V10	V11	V12	V13	V14	V15	V16	V17	V18	V19	V20	V21	V22	V23	V24	V25	V26	V27	V28	V29	V30
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	1	0	0	1	0	0	0	0	1	0	0
0	1	0	0	0	0	1	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	1	0	0	0
0	0	0	1	0	0	0	1	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	0
0	0	0	0	1	0	0	0	0	0	1	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0
0	0	0	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	1	0	0	0	0	0	0	0	1	0
0	0	0	1	0	0	0	0	0	1	0	0	0	0	0	0	0	1	0	0	0	1	0	0	0	0	0	0	0	0
0	0	0	0	0	1	0	0	1	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	1
0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0
0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	1	0	0	0	0	1	0	0	0	0	0	1	0
0	1	0	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	1	0	0	0	0	0	1	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0
1	0	0	0	0	0	0	0	0	1	0	0	1	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0
1	0	0	0	0	0	0	1	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0
0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	1	0	0	0	0	0	1
0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	1	0	0	0	0	1	0	0	0	0	0	0