diff --git a/R/reco_fun.R b/R/reco_fun.R
index 6edc728..27fa2a0 100755
--- a/R/reco_fun.R
+++ b/R/reco_fun.R
@@ -468,6 +468,19 @@ reco.proj_immutable <- function(base, cons_mat, cov_mat, immutable = NULL, ...){
     cli_abort("{.code max(immutable)} must be less or equal to {NCOL(base)}", call = NULL)
   }
 
+  # check immutable feasibility
+  # TODO: can proj_immutable2 be more stable than proj_immutable?
+  # Answer issue: https://github.com/danigiro/FoReco/issues/6#issue-2397642027 (@AngelPone)
+  cons_mat_red <- cons_mat[ , -immutable, drop = FALSE]
+  cons_vec <- apply(-cons_mat[ , immutable, drop = FALSE], 1, function(w)
+    rowSums(base[, immutable, drop = FALSE]%*%w))
+  if(length(immutable)>2){
+    check <- which(rowSums(cons_mat_red != 0) == 0)
+    if(length(check) > 0 && any(cons_vec[, check] > sqrt(.Machine$double.eps))){
+      cli_abort("There is no solution with this {.arg immutable} set.",  call = NULL)
+    }
+  }
+
   # Complete constraints matrix (immutable forecasts + linear constraints)
   imm_cons_mat <- .sparseDiagonal(NCOL(base))[immutable, , drop = FALSE]
   imm_cons_vec <- base[, immutable, drop = FALSE]
@@ -554,42 +567,6 @@ reco.strc_immutable <- function(base, strc_mat, cov_mat, immutable = NULL, ...){
     cli_abort("{.code max(immutable)} must be less or equal to {NCOL(base)}", call = NULL)
   }
 
-  # bts <- find_bts(strc_mat)
-  # id <- idrev <- 1:NCOL(cons_mat)
-  # if(any(!(immutable %in% bts))){
-  #   if(is.null(cons_mat)){
-  #     agg_mat <- strc_mat[-bts, , drop = FALSE]
-  #     cons_mat <- Matrix(0, nrow = NROW(agg_mat), ncol = NROW(strc_mat))
-  #     cons_mat[, bts] <- -agg_mat
-  #     cons_mat[, -bts] <- .sparseDiagonal(NROW(cons_mat))
-  #   }
-  #   id <- c(id[!(id%in%immutable) & !(id%in%bts)], id[!(id%in%immutable) & (id%in%bts)], immutable)
-  #   qr <- lcmat(cons_mat[,id, drop = FALSE], method = "rref", verbose = TRUE)
-  #   strc_mat <- rbind(qr$agg_mat, .sparseDiagonal(NCOL(qr$agg_mat)))
-  #   id <- id[qr$pivot]
-  #   idrev <- sort(id, index.return = TRUE)$ix
-  #   base <- base[, id , drop = FALSE]
-  #   cov_mat <- cov_mat[id, id , drop = FALSE]
-  #   bts <- find_bts(strc_mat)
-  #   immutable <- which(id %in% immutable)
-  #   if(any(!(immutable %in% bts))){
-  #     cli_abort("There is no solution with this {.arg immutable} set.", call = NULL)
-  #   }
-  # }
-  #
-  # S1 <- strc_mat[-immutable, -c(immutable-NROW(cons_mat)), drop = FALSE]
-  # S2 <- strc_mat[-bts, c(immutable-NROW(cons_mat)), drop = FALSE]
-  # S2u <- base[, immutable, drop = FALSE]%*%t(S2)
-  # base2 <- Matrix(base)
-  # base2[, !(1:NCOL(base2) %in% bts)] <- (base[, -bts] - S2u)
-  # reco_bts <- base2[, bts, drop = FALSE]
-  # cov_mat_red <- cov_mat[-immutable, -immutable, drop = FALSE]
-  # tmp <- reco.strc(base2[,-immutable, drop = FALSE], strc_mat = S1,
-  #                   cov_mat = cov_mat_red)[, (NROW(S1)-NCOL(S1)+1):NROW(S1), drop = FALSE]
-  # reco_bts[, !(bts %in% immutable)] <- tmp
-  # reco <- reco_bts %*% t(strc_mat)
-  # return(as.matrix(reco[,idrev, drop = FALSE]))
-
   # Code idea: https://github.com/AngelPone/chf
   bts <- find_bts(strc_mat)
   immutable <- sort(immutable)
@@ -608,7 +585,7 @@ reco.strc_immutable <- function(base, strc_mat, cov_mat, immutable = NULL, ...){
           i <- i - 1
           next
         }else{
-          cli_abort("Can not describe the hierarchy.",  call = NULL)
+          cli_abort("There is no solution with this {.arg immutable} set.",  call = NULL)
         }
       }
       determined <- determined[determined != immutable[i]]
@@ -633,7 +610,6 @@ reco.strc_immutable <- function(base, strc_mat, cov_mat, immutable = NULL, ...){
   return(as.matrix(reco))
 }
 
-
 transform_strc_mat <- function(strc_mat, bts){
   if (length(bts) != NCOL(strc_mat)){
     stop(simpleError(sprintf('length of basis set should be %d', NCOL(strc_mat))))
diff --git a/R/reco_opt.R b/R/reco_opt.R
index b2faefb..5c523e0 100755
--- a/R/reco_opt.R
+++ b/R/reco_opt.R
@@ -133,9 +133,15 @@ csrec <- function(base, agg_mat, cons_mat,
     if(NCOL(immutable) != 1){
       cli_abort("{.arg immutable} is not a vector.", call = NULL)
     }
-    if(max(immutable) > n){
+
+    if(max(immutable) >= n){
       cli_abort("{.code max(immutable)} must be less or equal to {n}", call = NULL)
     }
+
+    if(length(immutable) >= n){
+      # Answer issue: https://github.com/danigiro/FoReco/issues/6#issue-2397642027 (@AngelPone)
+      cli_abort("{.code length(immutable)} must be less than {n}", call = NULL)
+    }
   }
 
   # Compute covariance
@@ -309,7 +315,6 @@ terec <- function(base, agg_order, comb = "ols", res = NULL, tew = "sum",
 
   # Check immutable
   if(!is.null(immutable)){
-
     if(is.vector(immutable)){
       immutable <- matrix(immutable, ncol = length(immutable))
     }
@@ -326,8 +331,12 @@ terec <- function(base, agg_order, comb = "ols", res = NULL, tew = "sum",
                 do.call(c, sapply(m/kset, seq.int)) == x[2])
       }
     })
+
     if(is.null(immutable)){
       cli_warn("No immutable forecasts", call = NULL)
+    }else if(length(immutable) >= kt){
+      # Answer issue: https://github.com/danigiro/FoReco/issues/6#issue-2397642027 (@AngelPone)
+      cli_abort("The number of immutable constraints must be less than {kt}", call = NULL)
     }
   }
 
@@ -573,6 +582,12 @@ ctrec <- function(base, agg_mat, cons_mat, agg_order, comb = "ols", res = NULL,
       imm_mat[immutable[,1], col_id] <- 1
       imm_vec <- as(t(imm_mat), "sparseVector")
       immutable <- imm_vec@i
+
+      if(length(immutable) >= tmp$dim[["n"]]*tmp$dim[["kt"]]){
+        # Answer issue: https://github.com/danigiro/FoReco/issues/6#issue-2397642027 (@AngelPone)
+        cli_abort("The number of immutable constraints must be less than {tmp$dim[['n']]*tmp$dim[['kt']]}",
+                  call = NULL)
+      }
     }
   }