Skip to content

Commit bf206ad

Browse files
author
Ajay_kumar
committed
include separate flag for all the algorithms
1 parent 2d0b434 commit bf206ad

4 files changed

Lines changed: 233 additions & 50 deletions

File tree

‎src/Engine.cu‎

Lines changed: 33 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -40,14 +40,21 @@ Engine::Engine(DwnNetwork *myNetwork, ScenarioTree *myScenarioTree, SmpcConfigur
4040
cublasCreate(&handle);
4141
priceUncertaintyFlag = true;
4242
demandUncertaintyFlag = true;
43+
4344
if( algorithmName.compare("globalFbeAlgorithm") == 0){
44-
cout << "algorithm based on SMPC " << algorithmName <<endl;
4545
globalFbeFlag = true;
46-
46+
namaFlag = false;
47+
apgFlag = false;
48+
}else if(algorithmName.compare("namaAlgorithm") == 0){
49+
namaFlag = true;
50+
globalFbeFlag = false;
51+
apgFlag = false;
4752
}else{
48-
cout << "algorithm based on SMPC " << algorithmName <<endl;
53+
apgFlag = true;
54+
namaFlag = false;
4955
globalFbeFlag = false;
5056
}
57+
cout << "algorithm based on SMPC " << algorithmName <<endl;
5158

5259
_CUDA( cudaMalloc((void**)&devMatPhi, 2*nodes*nv*nx*sizeof(real_t)) );
5360
_CUDA( cudaMalloc((void**)&devMatPsi, nodes*nu*nv*sizeof(real_t)) );
@@ -142,13 +149,19 @@ Engine::Engine(SmpcConfiguration *smpcConfig){
142149
demandUncertaintyFlag = true;
143150

144151
if( algorithmName.compare("globalFbeAlgorithm") == 0){
145-
cout << " algorithm based on SMPC " << algorithmName <<endl;
146152
globalFbeFlag = true;
147-
153+
namaFlag = false;
154+
apgFlag = false;
155+
}else if(algorithmName.compare("namaAlgorithm") == 0){
156+
namaFlag = true;
157+
globalFbeFlag = false;
158+
apgFlag = false;
148159
}else{
149-
cout << " algorithm based on SMPC " << algorithmName <<endl;
160+
apgFlag = true;
161+
namaFlag = false;
150162
globalFbeFlag = false;
151163
}
164+
cout << "algorithm based on SMPC " << algorithmName <<endl;
152165

153166
_CUDA( cudaMalloc((void**)&devMatPhi, 2*nodes*nv*nx*sizeof(real_t)) );
154167
_CUDA( cudaMalloc((void**)&devMatPsi, nodes*nu*nv*sizeof(real_t)) );
@@ -1093,12 +1106,26 @@ bool Engine::getPriceUncertainty(){
10931106
bool Engine::getDemandUncertantiy(){
10941107
return demandUncertaintyFlag;
10951108
}
1109+
/*
1110+
* Flag to indicate status of the APG algorithm
1111+
*/
1112+
bool Engine::getApgFlag(){
1113+
return apgFlag;
1114+
}
1115+
10961116
/*
10971117
* Flag to indicate status of the globalFbe algorithm
10981118
*/
10991119
bool Engine::getGlobalFbeFlag(){
11001120
return globalFbeFlag;
11011121
}
1122+
/*
1123+
* Flag to indicate status of the NAMA algorithm
1124+
*/
1125+
bool Engine::getNamaFlag(){
1126+
return namaFlag;
1127+
}
1128+
11021129
/* SETTER'S IN THE ENGINE */
11031130
/*
11041131
* Option for uncertainty in price

‎src/Engine.cuh‎

Lines changed: 25 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -324,12 +324,25 @@ public:
324324
* status of the demand uncertanity
325325
*/
326326
bool getDemandUncertantiy();
327+
/*
328+
* bool to indicate APG algorithm
329+
* default 1 or true
330+
* globalFbe algorithm 1 or true
331+
*/
332+
bool getApgFlag();
327333
/*
328334
* bool to indicate gloablFbe algorithm
329-
* APG algorithm 0 or false - default
335+
* default 0 or false
330336
* globalFbe algorithm 1 or true
331337
*/
332338
bool getGlobalFbeFlag();
339+
/*
340+
* bool to indicate NAMA algorithm
341+
* default 0 or false
342+
* nama algorithm 1 or true
343+
*/
344+
bool getNamaFlag();
345+
333346
/* SETTER'S IN THE ENGINE */
334347
/*
335348
* Option for uncertainty in price
@@ -486,7 +499,7 @@ private:
486499
*/
487500
real_t *devSysUmin;
488501
/**
489-
* actuator/cotrol maximum
502+
* actuator/control maximum
490503
*/
491504
real_t *devSysUmax;
492505

@@ -642,10 +655,18 @@ private:
642655
*/
643656
bool demandUncertaintyFlag;
644657
/*
645-
* Flag to indicate the algorithm used (default 0 or false; uses APG
646-
* algorithm)
658+
* Flag to indicate the APG algorithm status (default 1 or true)
659+
*/
660+
bool apgFlag;
661+
/*
662+
* Flag to indicate the global FBE algorithm (default 0 or false)
647663
*/
648664
bool globalFbeFlag;
665+
/*
666+
* Flag to indicate the whether NAMA algorithm is used (default 0 or false; uses APG
667+
* algorithm)
668+
*/
669+
bool namaFlag;
649670
/**
650671
* cublas handler
651672
*/

‎src/SmpcController.cu‎

Lines changed: 121 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -79,6 +79,7 @@ SmpcController::SmpcController(string pathToConfigFile){
7979
factorStepFlag = false;
8080
simulatorFlag = true;
8181
bool globalFbeStatus = ptrMyEngine->getGlobalFbeFlag();
82+
bool namaStatus = ptrMyEngine->getNamaFlag();
8283

8384
vecPrimalInfs = new real_t[ptrMySmpcConfig->getMaxIterations()];
8485
vecValueFbe = new real_t[ptrMySmpcConfig->getMaxIterations()];
@@ -338,6 +339,8 @@ void SmpcController::initialiseAlgorithm(){
338339
_CUDA( cudaMemset(devVecPrevPsi, 0, nu*nodes*sizeof(real_t)) );
339340
_CUDA( cudaMemset(devVecGradientFbeXi, 0, 2*nx*nodes*sizeof(real_t)));
340341
_CUDA( cudaMemset(devVecGradientFbePsi, 0, nu*nodes*sizeof(real_t)));
342+
}else if (ptrMyEngine->getGlobalFbeFlag()){
343+
341344
}
342345
}
343346

@@ -891,13 +894,6 @@ void SmpcController::computeHessianOracalGlobalFbe(){
891894
}
892895
}
893896

894-
/*real_t* xDirHost = new real_t[nx*nodes];
895-
_CUDA( cudaMemcpy(xDirHost, devVecXdir, nx*nodes*sizeof(real_t), cudaMemcpyDeviceToHost) );
896-
897-
for( uint_t index = 0; index < 2*nx; index++)
898-
cout << xDirHost[index] << " ";
899-
900-
cout<< endl;*/
901897
// primalX = HxDir
902898
_CUBLAS(cublasSgemmBatched(ptrMyEngine->getCublasHandle(), CUBLAS_OP_N, CUBLAS_OP_N, 2*nx, 1, nx, &alpha, (const real_t**)
903899
ptrMyEngine->getPtrSysMatF(), 2*nx, (const real_t**)devPtrVecXdir, nx, &beta, devPtrVecPrimalXiDir, 2*nx, nodes) );
@@ -943,6 +939,85 @@ void SmpcController::computeGradientFbe(){
943939
// ptrMyEngine->getPtrSysMatG(), nu, (const real_t**)devPtrVecUDir, nu, &alpha, devPtrVecGradFbePsi, nu, nodes) );
944940
}
945941

942+
943+
/*
944+
* update the lbfgs buffer
945+
*/
946+
void SmpcController::updateLbfgsBuffer(){
947+
uint_t nodes = ptrMyEngine->getScenarioTree()->getNumNodes();
948+
uint_t nx = ptrMyEngine->getDwnNetwork()->getNumTanks();
949+
uint_t nu = ptrMyEngine->getDwnNetwork()->getNumControls();
950+
uint_t nDualXi = 2*nx*nodes;
951+
uint_t nDualPsi = nu*nodes;
952+
uint_t lbfgsBufferSize = ptrMySmpcConfig->getLbfgsBufferSize();
953+
real_t negAlpha = -1;
954+
real_t lbfgsBeta, lbfgsScale, invRhoVar;
955+
real_t normMatYk, normGrad, normMatSk;
956+
real_t gammaH;
957+
uint_t indexVecStart;
958+
real_t *devVecS;
959+
real_t *devVecY;
960+
961+
_CUDA( cudaMalloc((void**)&devVecS, (nDualXi + nDualPsi)*sizeof(real_t)) );
962+
_CUDA( cudaMalloc((void**)&devVecY, (nDualXi + nDualPsi)*sizeof(real_t)) );
963+
// vecS = y - yOld
964+
_CUDA( cudaMemcpy(devVecS, devVecXi, nDualXi*sizeof(real_t), cudaMemcpyDeviceToDevice) );
965+
_CUDA( cudaMemcpy(&devVecS[nDualXi], devVecPsi, nDualPsi*sizeof(real_t), cudaMemcpyDeviceToDevice) );
966+
_CUBLAS(cublasSaxpy_v2(ptrMyEngine->getCublasHandle(), nDualXi, &negAlpha, devVecPrevXi, 1, devVecS, 1));
967+
_CUBLAS(cublasSaxpy_v2(ptrMyEngine->getCublasHandle(), nDualPsi, &negAlpha, devVecPrevPsi, 1, &devVecS[nDualXi], 1));
968+
969+
if (ptrMyEngine->getGlobalFbeFlag()){
970+
// vecY = grad - gradOld
971+
_CUDA( cudaMemcpy(devVecY, devVecGradientFbeXi, nDualXi*sizeof(real_t), cudaMemcpyDeviceToDevice) );
972+
_CUDA( cudaMemcpy(&devVecY[nDualXi], devVecGradientFbePsi, nDualPsi*sizeof(real_t), cudaMemcpyDeviceToDevice) );
973+
_CUBLAS(cublasSaxpy_v2(ptrMyEngine->getCublasHandle(), nDualXi, &negAlpha, devVecPrevGradientFbeXi, 1, devVecY, 1));
974+
_CUBLAS(cublasSaxpy_v2(ptrMyEngine->getCublasHandle(), nDualPsi, &negAlpha, devVecPrevGradientFbePsi, 1, &devVecY[nDualXi], 1));
975+
// normGrad = norm(fbeGradient)
976+
_CUBLAS(cublasSnrm2_v2(ptrMyEngine->getCublasHandle(), nDualXi, devVecGradientFbeXi, 1, &lbfgsScale) );
977+
_CUBLAS(cublasSnrm2_v2(ptrMyEngine->getCublasHandle(), nDualPsi, devVecGradientFbePsi, 1, &normGrad) );
978+
normGrad = sqrt(lbfgsScale*lbfgsScale + normGrad*normGrad);
979+
}
980+
981+
if (ptrMyEngine->getNamaFlag()){
982+
983+
}
984+
// invRho = vecS'vecY
985+
_CUBLAS(cublasSdot_v2(ptrMyEngine->getCublasHandle(), nDualXi + nDualPsi, devVecS, 1, devVecY, 1, &invRhoVar) );
986+
// normMatYk = sqrt(vecY[:, lbfgsBufferCol]*vecY[:, lbfgsBufferCol])
987+
_CUBLAS(cublasSnrm2_v2(ptrMyEngine->getCublasHandle(), nDualXi + nDualPsi, devVecY, 1, &normMatYk) );
988+
// normMatSk = sqrt(matS[:, lbfgsBufferCol]*matS[:, lbfgsBufferCol])
989+
_CUBLAS(cublasSnrm2_v2(ptrMyEngine->getCublasHandle(), nDualXi + nDualPsi, devVecS, 1, &normMatSk) );
990+
991+
992+
if(normGrad < 1){
993+
normGrad = normGrad*normGrad*normGrad;
994+
}
995+
if( invRhoVar/(normMatSk * normMatSk) > 1e-6*normGrad){
996+
lbfgsBufferCol = 1 + (lbfgsBufferCol % lbfgsBufferSize);
997+
lbfgsBufferMemory = min(lbfgsBufferMemory + 1, lbfgsBufferSize);
998+
indexVecStart = lbfgsBufferCol*(nDualXi + nDualPsi);
999+
_CUDA( cudaMemcpy(&devLbfgsBufferMatY[indexVecStart], devVecY, (nDualXi + nDualPsi)*sizeof(real_t),
1000+
cudaMemcpyDeviceToDevice) );
1001+
_CUDA( cudaMemcpy(&devLbfgsBufferMatS[indexVecStart], devVecS, (nDualXi + nDualPsi)*sizeof(real_t),
1002+
cudaMemcpyDeviceToDevice) );
1003+
lbfgsBufferRho[lbfgsBufferCol] = 1/invRhoVar;
1004+
}else{
1005+
lbfgsSkipCount = lbfgsSkipCount + 1;
1006+
}
1007+
1008+
gammaH = invRhoVar/(normMatYk * normMatYk);
1009+
if( gammaH < 0 || abs(gammaH - lbfgsBufferHessian) == 0){
1010+
lbfgsBufferHessian = 1;
1011+
}else{
1012+
lbfgsBufferHessian = gammaH;
1013+
}
1014+
1015+
_CUDA( cudaFree(devVecS));
1016+
_CUDA( cudaFree(devVecY));
1017+
devVecS = NULL;
1018+
devVecY = NULL;
1019+
}
1020+
9461021
/*
9471022
* compute lbfgs direction
9481023
*/
@@ -1058,6 +1133,8 @@ void SmpcController::computeLbfgsDirection(){
10581133
devVecY = NULL;
10591134
}
10601135

1136+
1137+
10611138
/**
10621139
* function that compute the line search lbfgs update
10631140
*/
@@ -1221,7 +1298,7 @@ uint_t SmpcController::algorithmApg(){
12211298
}
12221299

12231300
/**
1224-
* This method executes the APG algorithm and returns the primal infeasibility.
1301+
* This method executes the global FBE algorithm and returns the primal infeasibility.
12251302
* @return primalInfeasibilty;
12261303
*/
12271304
uint_t SmpcController::algorithmGlobalFbe(){
@@ -1256,6 +1333,42 @@ uint_t SmpcController::algorithmGlobalFbe(){
12561333
return 1;
12571334
}
12581335

1336+
/**
1337+
* This method executes the NAMA algorithm and returns the primal infeasibilit
1338+
* @return primalInfeasibility;
1339+
**/
1340+
uint_t SmpcController::algorithmNama(){
1341+
1342+
initialiseAlgorithm();
1343+
initaliseLbfgBuffer();
1344+
1345+
uint_t maxIndex;
1346+
uint_t nx = ptrMyEngine->getDwnNetwork()->getNumTanks();
1347+
uint_t nu = ptrMyEngine->getDwnNetwork()->getNumControls();
1348+
uint_t nodes = ptrMyEngine->getScenarioTree()->getNumNodes();
1349+
1350+
for (uint_t iter = 0; iter < ptrMySmpcConfig->getMaxIterations(); iter++){
1351+
_CUDA( cudaMemcpy(devVecAcceleratedXi, devVecXi, 2*nx*nodes*sizeof(real_t), cudaMemcpyDeviceToDevice));
1352+
_CUDA( cudaMemcpy(devVecAcceleratedPsi, devVecPsi, nu*nodes*sizeof(real_t), cudaMemcpyDeviceToDevice));
1353+
solveStep();
1354+
proximalFunG();
1355+
computeFixedPointResidual();
1356+
computeGradientFbe();
1357+
if( iter > 0){
1358+
vecValueFbe[iter - 1] = computeValueFbe();
1359+
computeLbfgsDirection();
1360+
vecTau[iter - 1] = computeLineSearchLbfgsUpdate( vecValueFbe[iter - 1] );
1361+
}
1362+
dualUpdate();
1363+
_CUBLAS( cublasIsamax_v2(ptrMyEngine->getCublasHandle(), (2*nx + nu)*nodes, devVecResidual,
1364+
1, &maxIndex));
1365+
_CUDA( cudaMemcpy(&vecPrimalInfs[iter], &devVecResidual[maxIndex - 1], sizeof(real_t),
1366+
cudaMemcpyDeviceToHost));
1367+
}
1368+
1369+
return 1;
1370+
}
1371+
12591372
/**
12601373
* Invoke the SMPC controller on the current state of the network.
12611374
* This method invokes #updateStateControl, eliminateInputDistubanceCoupling

0 commit comments

Comments
 (0)