proc optmodel;
  number npts = 5000;
  number tz3 = 1.7;
  number pseed = 237.1531;
  number psize = 0.2;
  number pi = 4 * atan(1);
  number icr0 = 1 / npts;
  number incr = icr0 * 2 * pi;
  number xd{i in 1..npts} = ((1 + tz3^2) + cos(incr * (i - 1))) *
    cos(incr * (i - 1)) * (1 + psize * cos(incr * (i - 1) * pseed));
  number yd{i in 1..npts} = ((1 + tz3^2) + cos(incr * (i - 1))) *
    sin(incr * (i - 1)) * (1 + psize * cos(incr * (i - 1) * pseed));
  var z1 init 1.0;
  var z2 init 0.0;
  var z3 init 1.0;
  var x{i in 1..npts} init xd[i];
  var y{i in 1..npts} init yd[i];
  minimize f = sum {i in 1..npts} ( (x[i] - xd[i])^2  + (y[i] -
     yd[i])^2 );
  con cons1{i in 1..npts}: ((x[i] - z1)^2 + (y[i] - z2)^2)^2  - ((x[i]
     - z1)^2 + (y[i] - z2)^2) * (1 + z3^2)^2= 0.0;
   solve with ipnlp / tech=ipkrylov;
quit;