How can I simulate this toggle mechanism?

One way is to set up a DAE: See tutorial/DSolveExamplesOfDAEs and example/ModelConstrainedSystemsAsDAEs.

The constraint that the driver (bottom rotating link) has a fixed length is taken care of by initial conditions and the DE. There are two possible starting positions for the driven link. One might have to inspect the result of Solve to determine which is desired.

constraints = ComplexExpand@{
    (* SquaredEuclideanDistance[{a1[t], a2[t]}, {0, 0}] == r1^2, not needed *)
    SquaredEuclideanDistance[{b1[t], b2[t]}, {a1[t], a2[t]}] == rod^2,
    SquaredEuclideanDistance[{c1, c2}, {b1[t], b2[t]}] == r2^2};

eqn = {
      (* DE to rotate small link *)
    D[{a1[t], a2[t]}, t] == omega/r1 Cross[{a1[t], a2[t]}],
      (* initial conditions *)
    {a1[0], a2[0]} == {0, r1},
    {b1[0], b2[0]} == {b1[t], b2[t]} /.
      First @ Solve[constraints /. Thread[{a1[t], a2[t]} -> {0, r1}], {b1[t], b2[t]}]};

Example

Block[{omega = 1, r1 = 1, rod = 3, r2 = 2, c1 = 1, c2 = 3},
  sol = First @ NDSolve[{eqn, constraints},
     {a1[t], a2[t], b1[t], b2[t]},
     {t, 0, 2 Pi/omega}]];

With[{c1 = 1, c2 = 3},
 Manipulate[
  Graphics[
   {Line[{{0, 0}, {a1[t], a2[t]}, {b1[t], b2[t]}, {c1, c2}} /. sol /. t -> t0]},
   PlotRange -> {{-2, 6}, {-2, 6}}
   ],
  {t0, 0, Dynamic@sol[[1, 2, 0, 1, 1, 2]]}]
 ]

[The piston can be included via trigonometry or by adding another constraint.]

---

A random linkage:

Module[{sol, x, y, t},
 DynamicModule[{wheels = 6, omega = 1, joints, centers, constraints, 
   variables, initial, eqn, t0 = 0., direction = Forward},

  joints = 
   Sort@RandomReal[1, {wheels, 2}] + Table[{i, 0}, {i, wheels}];
  centers = 
   Sort@RandomReal[1, {wheels, 2}] + Table[{i, (-1)^i}, {i, wheels}];
  variables = 
   Transpose@{Array[x[#][t] &, wheels], Array[y[#][t] &, wheels]};
  constraints = ComplexExpand@{
     MapThread[
      SquaredEuclideanDistance[#1, #2] == 
        SquaredEuclideanDistance[#3, #4] &,
      {Most[variables], Rest[variables], Most[joints], Rest[joints]}
      ],
     MapThread[
      SquaredEuclideanDistance[#1, #3] == 
        SquaredEuclideanDistance[#2, #3] &,
      Rest /@ {variables, joints, centers}
      ]
     };
  initial = Thread /@ Thread[variables == joints /. t -> 0] // Flatten;
  eqn = {D[First@variables, t] == 
     omega/SquaredEuclideanDistance[First@joints, 
        First@centers] Cross[First@variables - First@centers]};

  Check[sol = 
    First@NDSolve[Flatten@{eqn, constraints, initial}, 
      Flatten[variables], {t, -2 Pi/omega, 2 Pi/omega}], 
   direction = ForwardBackward, {NDSolve::ndsz, NDSolve::ndcf}];

  Column[{
    Labeled[
     Manipulator[Dynamic[t0], First[x[1][t] /. sol /. t -> "Domain"], 
      AnimationDirection -> direction], "t", Left],
    Graphics[
     Dynamic@With[{pts = variables /. sol /. t -> t0},
       {{Red, Thin,
         MapThread[
          Circle, {centers, 
           MapThread[EuclideanDistance, {centers, joints}]}]},
        Line[pts], Line[Transpose[{centers, pts}]],
        EdgeForm[Black], Red, Disk[#, 0.1] & /@ pts
        }],
     PlotRange -> {{0, wheels + 2}, {-2, 2}}
     ]
    }]
  ]]

I did a solution with contour tracing on the distance function. It gets pretty unstable sometimes, but it's a fun question to experiment with interactivity.

DynamicModule[{p1 = {0, 2}, p2 = {1, 3}, angles = {0, 0}, distance, 
  grad, tangent}, 
 distance[a1_, a2_] := 
  Norm[{Cos@a1, Sin@a1} - (Norm[p2 - p1] {Cos@a2, Sin@a2} + p1)]; 
 grad = D[distance[a1, a2], {{a1, a2}}] /. Abs' -> Sign; 
 tangent = Cross@grad; 
 Column[{Dynamic[
    angles = 
     Mod[# /. FindRoot[distance @@ # == distance @@ angles, {t, 0}] &[
       grad t + angles + .05 tangent /. Thread[{a1, a2} -> angles]], 
      2 Pi]; Graphics[{{Dashed, Circle[], Circle[p1, Norm[p2 - p1]], 
       Locator@Dynamic@p1, Locator@Dynamic@p2}, Red, 
      Line[{{0, 0}, {Cos@#, Sin@#}, 
          p1 + Norm[p2 - p1] {Cos@#2, Sin@#2}, p1} & @@ angles]}, 
     PlotRange -> 5, Frame -> True]], 
   Dynamic@LocatorPane[Dynamic@angles, 
     ContourPlot[distance[a1, a2], {a1, 0, 2 Pi}, {a2, 0, 2 Pi}, 
      ContourLabels -> True]]}]]

You can do this without a NDSolve by calculating the distance from the follower cranks joint to the end of the driving cranks end. Then use this distance with law of cosines to calculate the deviation angle. This is also pretty easy to implement on ANY hardware capable of doing a ArcCos and Atan2 operation (note that in c atan2 parameters are swapped).

With[{p0 = {0, 0}, p3 = {1, 3}, r1 = 1, r2 = 2, r3 = 3},
 Module[{p1, p2, s, d, a1, a2},
  Animate[
   p1 = r1 {Cos[a], Sin[a]};
   s = p1 - p3;
   d = Norm[s]; 
   a1 = ArcTan[s[[1]], s[[2]]];
   a2 = a1 + ArcCos[( d^2 + r2^2 - r3^2)/(2 d r2)];
   p2 = p3 + r2 {Cos[a2], Sin[a2]};
   Graphics[
    Line[{p0, p1, p2, p3}],
    PlotRange -> {{-2, 6}, {-2, 6}}
    ],
   {a, 0, 2 Pi}]
  ]]

Results in the same thing as the accepted answer:

This solution is much much simpler, but won't work in 3D cases. The accepted answer will, and its a more flexible model. Finally, the position of the piston, which wasn't calculated by other posters as trivial is:

pp = {r1*Cos[a] + r4 + Sin[a]^2/(2 r4), 0};

And you get the following by fiddling a bit with the graphics elements: