{ optimización con restricciones }
unit miner2;
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;

//Acumula el valor de las restricciones violadas y se lo suma
//ponderado por rel al valor de la funcion objetivo
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);
		ProblemaPtr^[k].uv:= w;
		if w > 0 then m:= m + w;
	end;
	m:=	ProblemaPtr^[0].f(0, x) + m*rel;
	frel:= m;
end;

{Acumula los vetores gradientes de las restricciones que
esten violadas (ProblemaPtr^[k].uv > 0) y se lo suma ponderado
por rel al gradiente de la funcion objetivo.
grel(x) = g0(x) + rel * (g1(x) + g2(x) + ... + gNrs(x))
}
procedure grel(kf: integer; var grad, X: TVectR ); far;
var
	k: integer;
	mg: TVectR;
begin
	mg:= TVectR.Create_Init( X.N );
	mg.Ceros;
	for k:= 1 to NRs do
	begin
		if ProblemaPtr^[k].uv > 0 then
		begin
			ProblemaPtr^[k].g(k, grad, x);
			mg.SumRPV( 1{2*ProblemaPtr^[k].uv}, grad );
		end;
	end;
	ProblemaPtr^[0].g(0, grad, X);
	grad.SumRPV( rel, mg );
	mg.Free;
end;

//Obtiene memoria para el problema
procedure Init( NVariables, NRestricciones: integer );
begin
	NVs:= NVariables;
	NRs:= NRestricciones;
	GetMem( ProblemaPtr, SizeOF(TFichaFunc)*(NRs+1) );
end;

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

//Libera la memoria asociada al problema
procedure Done;
begin
	FreeMem( ProblemaPtr, SizeOF(TFichaFunc)*(NRs+1) );
end;

{Minimiza a traves del minero con control de paso en S la
funcion objetivo relajada ponderando el objetivo original
por relax
}
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 );
end;

{Llama al minimizar relajado con distintos valores para relax
paso y errs. relax aumenta con cada vuelta de la iteracion, por
lo que cada vez van a ir afectando mas las restricciones en el
calculo de frel y grel. paso disminuye con cada paso, con lo
que la busqueda se hace cada vez mas exhaustiva y errs
disminuye para permitir la busqueda con pasos mas chicos.
}
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:= TVectR.Create_Init(X.N);
	xa.copy(X);
	nitersrel:= 0;
	krel:= 0;
	dx:= pasoini;
	repeat
		inc(krel);
		inc(nitersrel);
		relax:= krel*sqrt(krel);
		pasorel:= {dx;} pasoini/sqrt(relax);
		if pasorel < errs then pasorel:= 1.1*errs;
		errsrel:= pasorel/10;
		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.Free;
end;

end.

