How can I minimize features of the trainded model?
There are several algorithms which can help you in a smart way.

Usually, those algorithms are used to tune the hyperparameters of a model, so this is what you will find in the tutorials/examples. In your case, you have to find a good set of features instead of a good set of hyperparameters, but the principle is the same.

My suggestions:

1) SMAC. This is based on Bayesian optimization. It's an iterative process where a proxy function is built and maximized:


the function to be optimized (your XGBoost model) is evaluated in a point (in the feature's hyperspace) where the optimizer believes it can find the maxima (or, in the very first iteration, in a point given by the user); 
the result is added to the set of all the evaluation points, and this set is used to build the proxy function;
the proxy function is maximized, and the coordinates of that maximum are believed to be the same where the original function will have a maximum too. 


Those three steps are repeated as much as you want. So, repeat from the first step;

It works both for continuous and for categorical features, and you can also impose some constraints between features.

Here an example for your case, in Python (code not tested):

from smac.configspace import ConfigurationSpace
from ConfigSpace.hyperparameters import UniformFloatHyperparameter, UniformIntegerHyperparameter
from smac.scenario.scenario import Scenario
from smac.facade.smac_facade import SMAC

#a continuous feature that you know has to lie in the [25 ~ 40] range
cont_feat = UniformFloatHyperparameter("a_cont_feature", 25., 40., default_value=35.)

#another continuous feature, [0.05 ~ 4] range
cont_feat2 = UniformFloatHyperparameter("another_cont_feature", 0.05, 4, default_value=1)


#a binary feature
bin_feat = UniformIntegerHyperparameter("a_bin_feature", 0, 1, default_value=1)

#the configuration space where to search for the maxima
cs = ConfigurationSpace()

cs.add_hyperparameters([cont_feat, cont_feat2, bin_feat])


# Scenario object
scenario = Scenario({"run_obj": "quality",   # we optimize quality
                     "runcount-limit": 1000,  # maximum function evaluations
                     "cs": cs,               # the configuration space
                     "cutoff_time": None
                     })

#here we include the XGBoost model
def f_to_opt(cfg):

    #here be careful! Your features need to be in the correct order for a correct evaluation of the XGB model
    features = {k : cfg[k] for k in cfg if cfg[k]}
    prediction = model.predict(features)

    return prediction


smac = SMAC(scenario=scenario, rng=np.random.RandomState(42),
        tae_runner=f_to_opt)
opt_feat_set = smac.optimize()

#the set of features which maximize the output
print (opt_feat_set)


2) dlib optimisation. This converge much faster than the previous. As disclaimer, I have to say that this is an algorithm which in principle works only with functions that fulfill a certain criteria, and XGBoost models as functions do not. But in the reality it turns out that this procedure works also for less stringent functions, at least in the cases I tried it. So maybe you want also give a try.

An example code:

import dlib

#here we include the XGBoost model. Note that we cannot use categorical/integer/binary features
def f_to_opt(cont_feat, cont_feat2):
    return model.predict([cont_feat, cont_feat2])


x,y = dlib.find_max_global(holder_table, 
                           [25, 0.05],  # Lower bound constraints on cont_feat and cont_feat2 respectively
                           [40, 4],    # Upper bound constraints on cont_feat and cont_feat2 respectively
                           1000)         # The number of times find_max_global() will call  f_to_opt