Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- nrow=16;
- ncol = 24;
- lambda = matrix(sample.int(100, size = ncol*nrow, replace = T),nrow,ncol);
- lambda = lambda - diag(lambda)*diag(x=1, nrow, ncol);
- y = rpois(ncol,lambda) + rtruncnorm(ncol,0,1,mean = 0, sd = 1);
- x = matrix (0, nrow, 1);
- x_A1 = y[1]+y[2]+y[3];
- x_A2 = y[4]+y[7]+y[3];
- x_B1 = y[4]+y[5]+y[6];
- x_B2 = y[11]+y[1];
- x_C1 = y[7]+y[8]+y[9];
- x_C2 = y[2]+y[5]+y[12];
- x_D1 = y[10]+y[11]+y[12];
- x_D2 = y[3]+y[6]+y[9];
- x_E1 = y[13]+y[14]+y[15];
- x_E2 = y[18]+y[19]+y[23];
- x_F1 = y[20]+y[21]+y[19];
- x_F2 = y[22]+y[16]+y[13];
- x_G1 = y[23]+y[22]+y[24];
- x_G2 = y[14]+y[17]+y[20];
- x_H1 = y[16]+y[17]+y[18];
- x_H2 = y[15]+y[21]+y[24];
- d <- c(x_A1, x_A2,x_B1, x_B2,x_C1, x_C2,x_D1, x_D2,x_E1,
- x_E2,x_F1, x_F2,x_G1, x_G2,x_H1, x_H2)
- x <- matrix(d, nrow, byrow=TRUE)
- A = matrix(c(1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, #x_A^1
- 0,0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0, #x_A^2
- 0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, #x_B^1
- 1,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0, #x_B^2
- 0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, #x_C^1
- 0,1,0,0,1,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0, #x_C^2
- 0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0, #x_D^1
- 0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, #x_D^2
- 0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0, #x_E^1
- 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,1,0, #x_E^2
- 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0, #x_F^1
- 0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,0,0,0,1,0,0, #x_F^2
- 0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0,0,0, #x_G^2
- 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1, #x_G^1
- 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0, #x_H^1
- 0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,1,0,0,1), #x_H^2
- nrow, ncol, byrow= TRUE)
- eval_f <- function( yhat ) {
- return( list( "objective" = norm((mean(yhat-y))^2, type = "2")))
- }
- # inequality constraint
- eval_g_ineq <- function( yhat ) {
- constr <- c(0 - yhat)
- return( list( "constraints"=constr ))
- }
- # equalities constraint
- eval_g_eq <- function( yhat ) {
- constr <- c( x-A%*%yhat )
- return( list( "constraints"=constr ))
- }
- x0 <- y
- #lower bound of control variable
- lb <- c(matrix (0, ncol, 1))
- local_opts <- list( "algorithm" = "NLOPT_LD_MMA",
- "xtol_rel" = 1.0e-7 )
- opts <- list( "algorithm" = "NLOPT_LD_AUGLAG",
- "xtol_rel" = 1.0e-7,
- "maxeval" = 1000,
- "local_opts" = local_opts )
- res <- nloptr( x0=x0,
- eval_f=eval_f,
- eval_grad_f = NULL,
- lb=lb,
- eval_g_ineq = eval_g_ineq,
- eval_g_eq=eval_g_eq,
- opts=opts)
- print(res)
- **#model <- list()
- #model$B <- A
- #model$obj <- norm((y-yhat)^2, type = "2")
- #model$modelsense <- "min"
- #model$rhs <- c(x,0)
- #model$sense <- c('=', '>=')
- #model$vtype <- 'C'
- #result <- gurobi(model, params)
- #print('Solution:')
- #print(result$objval)
- #print(result$yhat)**
Add Comment
Please, Sign In to add comment