from ortools.sat.python import cp_model
model = cp_model.CpModel()

n = 4

# Create s_i variables
S = {}
for i in range(n+1):
    S[i] = model.NewIntVar(0, n+1, f's_{i}')

# Reified constraints: eq[i, j] <-> s_i == j
eq = {}
for i in range(n+1):
    for j in range(n+1):
        eq[i, j] = model.NewBoolVar(f's_{i} == {j}')
        model.Add(S[i] == j).OnlyEnforceIf(eq[i, j])
        model.Add(S[i] != j).OnlyEnforceIf(eq[i, j].Not())
        
# s_i = number of occurrences of i in sequence
for i in range(n+1):
    model.Add(
        S[i] == sum(eq[j, i] for j in range(n+1))
    )
    
solver = cp_model.CpSolver()
if solver.Solve(model) == cp_model.OPTIMAL:
    print([f's_{i}={solver.Value(S[i])}' for i in range(n+1)])