Removing points from a triangular array without losing triangles

Python 3, n=8

import itertools
from ortools.sat.python import cp_model


def solve(n):
    model = cp_model.CpModel()
    solver = cp_model.CpSolver()
    cells = {
        (y, x): model.NewBoolVar(str((y, x)))
        for y in range(n) for x in range(y + 1)}
    triangles = [
            {cells[v] for v in ((y1, x1), (y2, x1), (y2, x1 + y2 - y1))}
            for (y1, x1) in cells.keys() for y2 in range(y1 + 1, n)]
    for t1, t2 in itertools.combinations(triangles, 2):
        model.AddBoolOr(t1.symmetric_difference(t2))
    model.Minimize(sum(cells.values()))
    solver.Solve(model)
    return len(cells) - round(solver.ObjectiveValue())


for n in itertools.count(2):
    print('a(%d) = %d' % (n, solve(n)))

Uses Google OR-Tools' CP-SAT solver.

After running for ~30 seconds, it outputs the following:

a(2) = 3
a(3) = 4
a(4) = 5
a(5) = 7
a(6) = 9
a(7) = 11
a(8) = 13

I didn't event try to wait for n=9 as it would probably take hours (the number of constraints grows like \$n^6\$). After less than 30 minutes of computation I found out that a(9)=15. I'm leaving my score as n=8 because at the moment the time constraints are unclear, but half an hour is probably too long.

How it works

Take two distinct equilateral triangles \$T_1\$ and \$T_2\$. To avoid ambiguity, there should be at least one bulb on a vertex belonging to exactly one of \$T_1\$ and \$T_2\$.

Thus the question may be rephrased as a SAT problem, with one constraint for every pair of triangles.

PS: I would very much like to include an example for n=8, but I'm having issues with the SAT solver which apparently wants to keep the solutions all for itself.

Getting the solutions from @Delfad0r's program

I extended @Delfad0r's program to output solutions. It also gives intermediate results, so you get output like this:

Solving n = 8:
a(8) >= 9
a(8) >= 10
a(8) >= 11
a(8) >= 12
a(8) >= 13
       o
      . o
     . o o
    . o o .
   o o . o o
  o o o o . .
 o . . o o o .
o . . o . o o o
a(8) = 13

Solving n = 9:
a(9) >= 10
a(9) >= 13
a(9) >= 14
a(9) >= 15
        o
       o o
      o . .
     o . o o
    . o . o o
   . o o o o o
  o o o . o . .
 o o o . . . o o
. o o o o o o . .
a(9) = 15

This computation took several hours.

If you get impatient and press Ctrl-C after some possibly non-optimal solution was found, the program will show that solution. So it doesn't take long to get this:

                   .
                  o o
                 . o o
                . o o o
               o o . o o
              o . o o o .
             o . o . o o o
            . o o o o o . o
           o . . o o o o o o
          o o o o o o o o o .
         o o . o o o o . o o o
        o o o o o o . o . o o o
       o . o o . o o o o o o o o
      o o o . o o o o o . o o o o
     o o o . o o o o o o o o . . o
    o o o o o o o o o o o . o . . o
   o o o o . o o o o . o o o o o . o
  o o o o o o o o . o o . . o o o o .
 o o o o . o o . o . o o o o o o . o o
o o . o o . o o o o . o o o . o o o o o
a(20) >= 42

Here is the extended program:

import itertools
from ortools.sat.python import cp_model

class ReportPrinter(cp_model.CpSolverSolutionCallback):

    def __init__(self, n, total):
        cp_model.CpSolverSolutionCallback.__init__(self)
        self.__n = n
        self.__total = total

    def on_solution_callback(self):
        print('a(%d) >= %d' %
              (self.__n, self.__total-self.ObjectiveValue()) )

def solve(n):
    model = cp_model.CpModel()
    solver = cp_model.CpSolver()
    cells = {
        (y, x): model.NewBoolVar(str((y, x)))
        for y in range(n) for x in range(y + 1)}
    triangles = [
            {cells[v] for v in ((y1, x1), (y2, x1), (y2, x1 + y2 - y1))}
            for (y1, x1) in cells.keys() for y2 in range(y1 + 1, n)]
    for t1, t2 in itertools.combinations(triangles, 2):
        model.AddBoolOr(t1.symmetric_difference(t2))
    model.Minimize(sum(cells.values()))
    print('Solving n = %d:' % n)
    status = solver.SolveWithSolutionCallback(model, ReportPrinter(n,len(cells)))
    if status == cp_model.OPTIMAL:
        rel = '='
    elif status == cp_model.FEASIBLE:
        rel = '>='
    else:
        print('No result for a(%d)\n' % n)
        return
    for y in range(n):
        print(' '*(n-y-1), end='')
        for x in range(y+1):
            print('.o'[solver.Value(cells[(y,x)])],end=' ')
        print()
    print('a(%d) %s %d' % (n, rel, len(cells) - solver.ObjectiveValue()))
    print()

for n in itertools.count(2):
    solve(n)