How to enforce specific elements in a vector to be in an optimization solution in Gekko

I'm trying to generate an optimal combination of p records in a vector of length n while simultaneously ensuring (constraining) that specific elements in the vector are included in the solution set p (based on a binary value for now)

I have the equation "simu_total_volume" below that does work for ensuring the solution set never exceeds p records, but I am unable to figure out how to modify this equation to ensure that specific elements i in the vector are included in the solution set (even if non-optimal). x7 is the binary vector that selects which elements are included based on p.

"Labor Day" is a vector of 0s except for one element that is equal to 1 (this corresponds to element I want to include in the solution set). I can enforce this vector to sum to 1 per the below, but am unsure how to integrate it into "simu_total_volume" so that the solution conforms to it.

Apologies for not including all relevant information for reproducibility, but the full solution is very large.

simu_total_volume = [m.Intermediate((
(m.max2(0,base_volume[i]*(m.exp(total_vol_fedi[i])-1)) * x3[i] +
m.max2(0,base_volume[i]*(m.exp(total_vol_feao[i])-1)) * x4[i] +
m.max2(0,base_volume[i]*(m.exp(total_vol_diso[i])-1)) * x5[i] +
m.max2(0,base_volume[i]*(m.exp(total_vol_tpro[i])-1)) * x6[i]) + base_volume[i]) * x7[i]) for i in range(n)]

labor_day = [m.Intermediate(x8[i] * el_cppg['holiday_labor_day_flag'].values[i]) for i in range(n)]

#Require labor day to be in output
m.Equation(sum(labor_day) == 1)

#Limit max output
m.Equation(sum(x7)<=p)

m.Maximize(m.sum(simu_total_volume))

m.options.SOLVER=1

try:
    m.solve(disp = True)
except:
    continue

There are multiple ways to constrain particular elements within an array. Here is a complete example that optimizes the elements of X to minimize the total cost. Each element X can be [0,1] and two can be selected with sum(X)==2.

from gekko import GEKKO
m = GEKKO(remote=False)
X = m.Array(m.Var,7,lb=0,ub=1,integer=True)
c = [1.2,0.95,1.3,1.0,0.8,1.25,1.4]
m.Equation(sum(X)==2)
m.Minimize(sum([X[i]*c[i] for i in range(7)]))
m.options.SOLVER=1
m.solve()
print(f'X: {X}')

The solver picks the two lowest that correspond to c elements of 0.95 and 0.8:

c = [1.2,0.95,1.3,1.0,0.8,1.25,1.4]
X: [[0.0] [1.0] [0.0] [0.0] [1.0] [0.0] [0.0]]

Here are a few ways to constrain the solution such as enforcing that the last element is always selected:

Add an equation

m.Equation(X[-1]==1)

Set the upper and lower bounds to the specified solution

X[-1].lower=1
X[-1].upper=1

Use the m.fix() function

m.fix(X[-1],1)

Add an objective as a soft constraint

Use this method if adding a hard constraint gives an infeasible solution. This will encourage the selection of a preferred option, but won't enforce it if the equations aren't satisfied.

m.Minimize(100*(X[-1]-1)**2)

Results

All of these methods return the correct solution that selects the last element (not optimal) and the least costly element.

X: [[0.0] [0.0] [0.0] [0.0] [1.0] [0.0] [1.0]]