With some diffidence (because there appears to be a Mathematica bug: see below), I would like to offer an answer in the spirit of the OP's original attempt to solve the problem algebraically.
Solution
This problem can be formulated as a binary integer linear program. The reformulation represents the square (or more generally, a rectangle as implemented below) as a vector. Its coefficients represent occupancy: $0$ for empty, $1$ for full. In this representation, every possible tile placement corresponds to a unique vector indicating the cells occupied by the tile.
Let's begin by generating the vectors corresponding to all possible placements of tiles in all possible orientations. It will be convenient to specify as little as possible, beginning with a set of unique tile shapes.
As a running example I use the L, C, T, and Z tiles:
tileSet = {{{1, 1, 1}, {1, 0, 0}}, {{1, 1, 1}, {1, 0, 1}}, {{1, 1, 1}, {0, 1, 0}}, {{0, 1, 1}, {1, 1, 0}}};
ArrayPlot[#, Frame -> False, ImageSize -> 10 Max[Length[First[#]]]] & /@ tileSet
(There is no restriction that tiles have a common bounding box, but--to make sure they can be properly positioned--it is important that their matrix descriptions not have any zero padding along outside rows or columns.)
We need to rotate (and, if we wish, flip) all these tiles. An elegant way (perhaps not the simplest) is to generate the group action abstractly in terms of generators (alpha and beta).
apply[s_List, alpha] := Reverse /@ s; (* Horizontal reflection *)
apply[s_List, beta] := Transpose[s]; (* Diagonal reflection *)
apply[s_List, g_List] := Fold[apply, s, g]; (* Group composition *)
group = FoldList[Append, {}, Riffle[ConstantArray[alpha, 4], beta]];
(Riffle and ConstantArray are just a trick, specific to Coxeter groups, to generate a unique set of group elements. If you don't want to allow reflections, then eliminate the four orientation-reversing group elements, which are the ones using odd numbers of generators.)
Apply the group action to the tiles and eliminate any duplicates:
tiles = Union[Flatten[Outer[apply[#1, #2] &, tileSet, group, 1], 1]];
ArrayPlot[#, Frame -> False, ImageSize -> 10 Max[Length[First[#]]]] & /@ tiles
These tiles need to be converted from arrays (suitable for graphical display) into vectors (for input to the linear program). Because these vectors represent all possible positions (translations) of the oriented tiles into the square, here is where we specify the square's dimensions, $m$ rows by $n$ columns. I use $6$ by $9$ for this example (which is a significant size). makeAllTiles creates all possible positions of a single tile within a rectangle; mapping it over the set of all oriented tiles and flattening the result gives what we want. (No need any more to remove duplicates, because there is no possibility of an accidental match.)
makeAllTiles[tile_, m_Integer, n_Integer] :=
With[{m0 = Length[tile], n0 = Length[First[tile]]},
Flatten[Table[ArrayPad[tile, {{i, m - m0 - i}, {j, n - n0 - j}}],
{i, 0, m - m0}, {j, 0, n - n0}], 1]];
allTiles = Flatten[makeAllTiles[#, 6, 9] & /@ tiles, 1];
ArrayPlot[options = Transpose[Flatten /@ allTiles], ImageSize -> 800]
Each possible position of each possible orientation of each tile shape is represented by a column vector in this array. The rows correspond to cells in the rectangle.
It's now surprisingly straightforward to complete the setup: we seek a subset of columns summing to the vector representing a filled square: all the coefficients in this sum must be $1$'s. The subset can be picked out by an indicator vector $x$. We therefore require its entries to be integral, bounded below by $0$, and bounded above by $1$. A feasible vector, let's say, will be one where the corresponding sum fills the square, but sometimes fills its cells multiple times (due to overlap of tiles): that is, a lower bound $b$ for the result should be a $1$ in each cell of the square. Finally, we minimize the sum of the cell values: if a solution exists, this sum will equal the number of squares. The sum is a linear combination with unit coefficients, represented as a vector $c$.
b = ConstantArray[1, Length[options]]; (* Lower bounds on cell values *)
c = ConstantArray[1, Length[First[options]]]; (* Sum of cell values *)
lu = ConstantArray[{0, 1}, Length[First[options]]]; (* 0-1 constraints *)
x = Quiet[LinearProgramming[c, options, b, lu, Integers, Tolerance -> 0.0005],
{LinearProgramming::lpip, LinearProgramming::lpsnf}];
If[!ListQ[x] || Max[options.x] > 1, x = {}];
solution = allTiles[[Select[x Range[Length[x]], # > 0 &]]];
ArrayPlot[z = Sum[i solution[[i]], {i, 1, Length[solution]}],
PlotRange -> {0, Max[z]}, Mesh -> True, ColorFunction -> "Rainbow"]
Solution is a list of tiles (in their orientations and positions).
Limitations
LinearProgramming (apparently) can return only one solution even when many are optimal. To explore all solutions, or a range of solutions, may require different methods. (There are still some tricks we can play: e.g., by permuting the order of columns in Options, we might wind up with a different solution.)
As with the systematic searching methods, when this problem gets big, the computing time grows rapidly. You can see this is so, just by considering the quadratic growth in the list of possible positions of the tiles as the square's dimensions increase. Computing time for the example is around one second (increasing as Tolerance decreases).
Finally, I had to fudge something to get this solution: LinearProgramming, with its default setting for Tolerance, returned a non-optimal solution! (options.x had some 2's in it, indicating some double overlaps.) By starting with Tolerance set to $0.5$ and gradually decreasing it, within a few steps a valid solution was obtained. (With the value set to 10^-1, though, I raised an Assert in the underlying C code :-(. That is a bug, at least in version 8.0.0.0.)
