# objective function
eval_f0 <- function(x0, x1, x2, x3, x4, x5){
x1 <- matrix(x0[1:nrow(trd)], ncol = 1)
return( t(x1)%*%x1)
}
# constraint function
eval_inque <- function(x0, x1, x2, x3, x4, x5) {
x2 <- matrix(x0[nrow(trd)+1:nrow(trd)*2], ncol = 1)
x3 <- matrix(x0[nrow(trd)*2+1:nrow(trd)*2+24], ncol = 1)
x4 <- matrix(x0[nrow(trd)*2+25:nrow(trd)*3+25], ncol = 1)
x5 <- matrix(x0[nrow(trd)*3+25:nrow(trd)*4+25], ncol = 1)
return(rbind(c(-1*x2),
c((trd%*%x3-x4%*%Price+x1)-x5),
c(-1*x5),
c(x2-1)))
}
# constraint function
eval_que <- function(x0, x1, x2, x3, x4, x5) {
x2 <- matrix(x0[nrow(trd)+1:nrow(trd)*2], ncol = 1)
x3 <- matrix(x0[nrow(trd)*2+1:nrow(trd)*2+24], ncol = 1)
x4 <- matrix(x0[nrow(trd)*2+25:nrow(trd)*3+25], ncol = 1)
x5 <- matrix(x0[nrow(trd)*3+25:nrow(trd)*4+25], ncol = 1)
return(rbind(c(x2-Market),
c(t(x2%*%(x5-(trd%*%x3-x4%*%Price+x1)))),
c(t(x5)%*%(1-x2))))
}
res1 <- nloptr( x0=rep(1, nrow(trd)*4+nrow(trd)),
eval_f = eval_f0,
eval_g_ineq = eval_inque,
eval_g_eq = eval_que,
lb = rep(0, nrow(trd)*4+nrow(trd)),
ub = rep(Inf, nrow(trd)*4+nrow(trd)),
opts = list("algorithm"="NLOPT_LD_SLSQP"),
x1 = matrix(nrow(trd), ncol = 1),
x2 = matrix(nrow(trd), ncol = 1),
x3 = matrix(23, ncol = 1),
x4 = matrix(nrow(trd), ncol = 1),
x5 = matrix(nrow(trd), ncol = 1)
)
IyBvYmplY3RpdmUgZnVuY3Rpb24KZXZhbF9mMCA8LSBmdW5jdGlvbih4MCwgeDEsIHgyLCB4MywgeDQsIHg1KXsKICB4MSA8LSBtYXRyaXgoeDBbMTpucm93KHRyZCldLCBuY29sID0gMSkKICByZXR1cm4oIHQoeDEpJSoleDEpCn0KIyBjb25zdHJhaW50IGZ1bmN0aW9uCmV2YWxfaW5xdWUgPC0gZnVuY3Rpb24oeDAsIHgxLCB4MiwgeDMsIHg0LCB4NSkgewogIHgyIDwtIG1hdHJpeCh4MFtucm93KHRyZCkrMTpucm93KHRyZCkqMl0sIG5jb2wgPSAxKQogIHgzIDwtIG1hdHJpeCh4MFtucm93KHRyZCkqMisxOm5yb3codHJkKSoyKzI0XSwgbmNvbCA9IDEpCiAgeDQgPC0gbWF0cml4KHgwW25yb3codHJkKSoyKzI1Om5yb3codHJkKSozKzI1XSwgbmNvbCA9IDEpCiAgeDUgPC0gbWF0cml4KHgwW25yb3codHJkKSozKzI1Om5yb3codHJkKSo0KzI1XSwgbmNvbCA9IDEpCiAgcmV0dXJuKHJiaW5kKGMoLTEqeDIpLAogICAgICAgICAgICAgICBjKCh0cmQlKiV4My14NCUqJVByaWNlK3gxKS14NSksCiAgICAgICAgICAgICAgIGMoLTEqeDUpLAogICAgICAgICAgICAgICBjKHgyLTEpKSkKfQojIGNvbnN0cmFpbnQgZnVuY3Rpb24KZXZhbF9xdWUgPC0gZnVuY3Rpb24oeDAsIHgxLCB4MiwgeDMsIHg0LCB4NSkgewogIHgyIDwtIG1hdHJpeCh4MFtucm93KHRyZCkrMTpucm93KHRyZCkqMl0sIG5jb2wgPSAxKQogIHgzIDwtIG1hdHJpeCh4MFtucm93KHRyZCkqMisxOm5yb3codHJkKSoyKzI0XSwgbmNvbCA9IDEpCiAgeDQgPC0gbWF0cml4KHgwW25yb3codHJkKSoyKzI1Om5yb3codHJkKSozKzI1XSwgbmNvbCA9IDEpCiAgeDUgPC0gbWF0cml4KHgwW25yb3codHJkKSozKzI1Om5yb3codHJkKSo0KzI1XSwgbmNvbCA9IDEpCiAgcmV0dXJuKHJiaW5kKGMoeDItTWFya2V0KSwKICAgICAgICAgICAgICAgYyh0KHgyJSolKHg1LSh0cmQlKiV4My14NCUqJVByaWNlK3gxKSkpKSwKICAgICAgICAgICAgICAgYyh0KHg1KSUqJSgxLXgyKSkpKQp9CgoKcmVzMSA8LSBubG9wdHIoIHgwPXJlcCgxLCBucm93KHRyZCkqNCtucm93KHRyZCkpLAogICAgICAgICAgICAgICAgZXZhbF9mID0gZXZhbF9mMCwKICAgICAgICAgICAgICAgIGV2YWxfZ19pbmVxID0gZXZhbF9pbnF1ZSwKICAgICAgICAgICAgICAgIGV2YWxfZ19lcSA9IGV2YWxfcXVlLAogICAgICAgICAgICAgICAgbGIgPSByZXAoMCwgbnJvdyh0cmQpKjQrbnJvdyh0cmQpKSwKICAgICAgICAgICAgICAgIHViID0gcmVwKEluZiwgbnJvdyh0cmQpKjQrbnJvdyh0cmQpKSwKICAgICAgICAgICAgICAgIG9wdHMgPSBsaXN0KCJhbGdvcml0aG0iPSJOTE9QVF9MRF9TTFNRUCIpLAogICAgICAgICAgICAgICAgeDEgPSBtYXRyaXgobnJvdyh0cmQpLCBuY29sID0gMSksCiAgICAgICAgICAgICAgICB4MiA9IG1hdHJpeChucm93KHRyZCksIG5jb2wgPSAxKSwKICAgICAgICAgICAgICAgIHgzID0gbWF0cml4KDIzLCBuY29sID0gMSksCiAgICAgICAgICAgICAgICB4NCA9IG1hdHJpeChucm93KHRyZCksIG5jb2wgPSAxKSwKICAgICAgICAgICAgICAgIHg1ID0gbWF0cml4KG5yb3codHJkKSwgbmNvbCA9IDEpCiAgICAgICAgICAgICAgICAp