{
Tecnica INEXTOR
===============


optimización con restricciones

las restricciones se modifican a los efectos de hacer que
si est n activas den algo de costo que tiende a cero al alejarse
de la restriccion y que es exponencialmente creciente aumentar la
violacion.

Se supone que las restrcciones estan planteadas de manera que el valor
de la funcion es representativo de la distancia del punto a la forntera.

 La restricci˘n modificada es:

 1/relax * exp( relax * r(x))


 El gradiente de la restricci˘n modificada es:

 exp(relax* r(x)) *  grad(r(X))



}
unit miner3;
interface
uses
	xMatDefs, miner, MatReal;

procedure Init( NVariables, NRestricciones: integer );
procedure DefFunc(
	irestric: integer;
	frestric: TFunc_fx;    { dominio= frestric(X) < 0 }
	grestric: TGrad_fx );

procedure MinRelax(
	relax: NReal;
	var X: TVectR;
	paso: NReal;
	errs: NReal;
	MaxNIters: longint;
	var NIters: longint;
	var ValCosto: NReal;
	var Convergio: boolean );

procedure Minimizar(
	var X: TVectR;
	pasoini: NReal;
	errs: NReal;
	MaxNItersRel: longint;
	var NItersRel: longint;
	MaxNIters: longint;
	var NIters: longint;
	var ValCosto: NReal;
	var Convergio: boolean );

procedure Done;

implementation

const
	MaxNRestric= 1000;

type
	TFichaFunc= record
		uv: NReal;
		f: TFunc_fx;
		g: TGrad_fx;
	end;

  PProblema= ^TProblema;
  TProblema= array[0..MaxNRestric] of TFichaFunc;

var
	ProblemaPtr: PProblema;
	rel: NReal;
	NVs, NRs: integer;
	FiltCoefB: NReal;

//retorna exp(x) limitando x a 11000 o a 80 segun este
//definido N+ o no
function expsat( x: NReal ): NReal;
begin
	{$IFOPT N+}
	if x > 11000 then x:= 11000;
	{$ELSE}
	if x > 80 then x:= 80;
	{$ENDIF}
	expsat:= exp(x);
end;

{Acumula los valores de las restricciones r(x) como
exp(r(x))/rel y se los suma al valor de la funcion objetivo
ponderado por rel.
En uv guarda solamente exp(rel*r(x)) porque es lo que se necesita
para calcular el gradiente
}
function frel(kf: integer; var X: TVectR ): NReal; far;
var
	k: integer;
	m, w: NReal;
begin
	m:=0;
	for k:= 1 to NRs do
	begin
		w:= ProblemaPtr^[k].f(k, x);
		w:= expsat(rel*w);
		ProblemaPtr^[k].uv:= w;
		m:= m + w/rel;
	end;
	m:=	ProblemaPtr^[0].f(0, x) + m*rel;
	frel:= m;
end;

{Acumula los gradientes de las restricciones modificadas como
exp(rel*r(x))*grad(r(x)) y se lo suma al gradiente de la funcion objetivo
ponderado por rel
}
procedure grel(kf: integer; var grad, X: TVectR ); far;
var
	k: integer;
	mg: TVectR;
begin
	mg.init( X.N );
	mg.Ceros;
	for k:= 1 to NRs do
	begin
		ProblemaPtr^[k].g(k, grad, x);
		mg.SumRPV( ProblemaPtr^[k].uv, grad );
	end;
	ProblemaPtr^[0].g(0, grad, X);
	grad.SumRPV( rel, mg );
	mg.done;
end;

procedure Init( NVariables, NRestricciones: integer );
begin
	NVs:= NVariables;
	NRs:= NRestricciones;
	GetMem( ProblemaPtr, SizeOF(TFichaFunc)*(NRs+1) );
end;

procedure DefFunc(
	irestric: integer;
	frestric: TFunc_fx;
	grestric: TGrad_fx );
begin
	ProblemaPtr^[irestric].f:= frestric;
	ProblemaPtr^[irestric].g:= grestric;
	ProblemaPtr^[irestric].uv:= 0;
end;

procedure Done;
begin
	FreeMem( ProblemaPtr, SizeOF(TFichaFunc)*(NRs+1) );
end;

procedure MinRelax(
	relax: NReal;
	var X: TVectR;
	paso: NReal;
	errs: NReal;
	MaxNIters: longint;
	var NIters: longint;
	var ValCosto: NReal;
	var Convergio: boolean );

begin
	rel:= relax;
{
	MinerCPS(	frel, grel, X, paso, errs, MaxNITers, NIters, ValCosto, Convergio );
}
{
	MinerCPX(	frel, grel, X, paso, errs, MaxNITers, NIters, ValCosto, Convergio );
}

	MinerFPP(	frel, grel, X, paso, errs, MaxNITers, NIters, ValCosto, Convergio, FiltCoefB );

end;



procedure Minimizar(
	var X: TVectR;
	pasoini: NReal;
	errs: NReal;
	MaxNItersRel: longint;
	var NItersRel: longint;
	MaxNIters: longint;
	var NIters: longint;
	var ValCosto: NReal;
	var Convergio: boolean );
var
	krel: integer;
	relax: NReal;
	xa: TVectR;
	errsrel, pasorel: NReal;
	dx: NReal;
begin
	xa.init(X.N);
	xa.copy(x);
	nitersrel:= 0;
	krel:= 0;
	dx:= pasoini;
	repeat
		inc(krel );
		inc(nitersrel);
		relax:= krel*sqrt(krel);
		pasorel:= {dx;} pasoini/relax;
		if pasorel < errs then pasorel:= 1.1*errs;
		errsrel:= pasorel/20;
		if errsrel< errs then errsrel:= errs;
		MinRelax(
				relax,
				X,
				pasorel,
				errsrel,
				MaxNIters,
				NIters,
				ValCosto,
				Convergio);
		Xa.SumRPV( -1, X);
		dx:= Xa.NormEuclid;
		Xa.copy( X );

	write('Rx: ', relax:12,'  NI: ', NIters:4);
	if convergio then
		writeln( '  x1: ',X.e(1):8:-4, '  x2: ', X.e(2):8:-4,' f(x): ', ValCosto:8:-4 )
	else
		writeln( '  NO CONVERGIO ');
	until (dx < errs )or( nitersrel > MaxNItersRel ) or not Convergio;
	xa.done;
end;

begin
	FiltCoefB:= 0.1;
end.

