Some additional explainations:

In my question, the values of n and lst are not arbitrary integers. They have the following relationships and restrictions:

1. They are all positive integers.
2. n>(#[[1]]&/@lst//Max)
3. The list elements are generated by Phi(x), where Phi is [Euler's totient function](https://en.wikipedia.org/wiki/Euler%27s_totient_function), and x is a prime or power of prime with the condition x - 1 <= n.

See following confirmation:

In[26]:= n>(#[[1]]&/@lst//Max)
(*
Here, #[[2]] must a prime:
*)
ForAll[#[[1]]==EulerPhi[#[[2]]^#[[3]]]&/@lst,True]

Out[26]= True

Out[27]= True

I want to know, under the constraints of the above conditions, whether this problem can be NP-Complete or even simpler?

Regards, Zhao

Thank you for your actual data. For that, at most 20 subsets can be in any solution. Perhaps your original code with

sub=Subsets[pk,{2,20}]

might be interesting to measure. That should be significantly faster because of trying far far fewer subsets. Please let me know how much that helps if you take the time to test that.

I am sorry that my idea did not provide sufficient speed-up for you.

Since you are using extremely long lists and you only want subsets where the first elements sum to n then it may be possible that you do not need to look at almost all possible subsets.

I sorted your list of lists on the first element. I created a recursive function to generate subsets. That will be much slower than the internal Mathematica Subsets function. But before each recursive step I looked at the total thus far and did not recurse if the total would be greater than n. That meant I could skip all remaining sets in the sorted list. For your sample data set that meant I only needed to look at 106 of the possible 1024 cases. This method was approximately 4.6 times faster than your original code for that sample data. Someone else might be able to see how to make method this even faster.

For much larger lists this might avoid an even greater amount of work. I do not know how easily this might be done with the code using IntegerDigits.

This example only generates the subset lists and does not impose your second condition.

Please test this very carefully before you depend on it.

AbsoluteTiming[pk=SortBy[{{2,3,1},{6,3,2},{18,3,3},{4,5,1},{20,5,2},{6,7,1},{10,11,1},{12,13,1},{16,17,1},{18,19,1}},First];
n=22;
count=0;
subsets[subset_,total_,list_]:=(
  Print[count++];
  If[total==n,Sow[subset]];
  If[list!={}&&total+list[[1,1]]<=n,
    subsets[subset,total,Rest[list]];
    subsets[Join[subset,{list[[1]]}],total+list[[1,1]],Rest[list]]
  ]);
Reap[subsets[{},0,pk]][[2,1]]]

I didn't analyze your algorithm, but I don't think it finds all of the possibilities. For example, your output is missing {{2, 3, 1}, {4, 5, 1}, {6, 3, 2}, {10, 11, 1}}.

I came up with an algorithm that is faster and that finds more results. It uses IntegerPartitions instead of Subsets. However it typically takes over a minute when length of pk is 30. How long will your pk be? Do you need to find absolutely all solutions?